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ABSTRACT 


We develop models and solution methodologies to solve the discrete-time path- 
optimization problem where a single or multiple searchers look for a moving target 
in a hnite set of cells. The single searcher is constrained by maximum limits on 
the consumption of several resources such as time, fuel, and risk along any path. 
We develop a specialized branch-an-bound algorithm for this problem that utilizes 
several new network reduction procedures as well as new bounding technique based 
on Lagrangian relaxation and network expansion. The resulting algorithm is quite 
efficient and promising. For the multiple searchers, an optimal set of paths (search 
plan) is determined by taking advantage of the cooperative search effect. We present 
a new exact algorithm and two new heuristics to hnd an optimal or near-optimal 
search plan. One of the heuristics is based on the cross-entropy method and is found 
to perform well for a broad range of problem instances. The exact algorithm deals with 
the specihc case of homogeneous searchers and is based on outer approximations by 
several new cutting planes. In addition, we prove that under certain assumptions the 
path-optimization problem becomes equivalent to a large-scale linear mixed-integer 
program. 
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EXECUTIVE SUMMARY 


This dissertation develops models and solution methodologies to solve the dis¬ 
crete time-and-space path-optimization problems where a single searcher or multiple 
searchers look for a non-evading moving target. We extended the single searcher 
problem (SSP), generally called the optimal searcher path problem, to realistic situ¬ 
ation where: (1) the searcher needs to change its flight altitude during a mission to 
balance sensor quality and risk from ground threats; (2) the searcher is subject to 
constraints on multiple resources such as fuel consumption and risk exposure along 
any path; and (3) the search effect depends on the current and previous locations of 
the searcher. The latter situation may occur if moving from some search sector to 
a new sector results in a lower search effect due to ineffective search during transit. 
We refer to this extended SSP as the resource-constrained single searcher problem 
(RSSP). 

We present new bounding techniques (static bound and directional static 
bound) and network reduction procedure in a branch-and-bound algorithm to opti¬ 
mally solve the single searcher problem (SSP). The static bound simplihes the bound 
calculation substantially at the expense of a weaker bound. The directional static 
bound improves the weaker static bound. The network reduction procedure takes 
advantage of the initial positions of searcher and target, and eliminates parts of the 
branch-and-bound tree while retaining a optimal solution. The resulting algorithm 
solves a problem instance with 15 by 15 cells and a time horizon of 20 optimally in 
less than 4 seconds on average, which is a signihcant reduction from the 9 minutes 
required by a state-of-the-art algorithm. 

We extend our branch-and-bound algorithm for the SSP to treat the RSSP 
and construct a new bounding technique based on Lagrangian relaxation. In addition, 
we develop novel network reduction techniques using dominance tests. The resulting 
algorithm solves a realistic problem instance (10 by 10 cells, two altitudes, time 
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horizon 40, and two resource constraints), on average, in less than 10 minutes. 

We extend the specialized branch-and-bound algorithm to the case of multi¬ 
ple searchers and also develop two new heuristics (static bound heuristic and cross¬ 
entropy heuristic). In the context of search problems, this dissertation appears to be 
the first one which utilizes the cross-entropy method. Among these three algorithms, 
the cross-entropy heuristic performs best for a broad range of problem instances. 
While the branch-and-bound algorithm and the static bound heuristic are limited 
to small problem instances, the cross-entropy heuristic obtains good approximate 
solutions of moderately sized instances (15 by 15 cells, time horizon 18, and three 
searchers) in about 10 minutes. 

For the case with multiple identical searchers, we construct an exact algorithm 
based on outer approximations by cutting planes. We introduce several new cutting 
planes (multiple-cut, strong-cut, relative-cut, and symmetric-cut), combine them, and 
develop a specialized outer approximation algorithm. For a moderate-size problem 
instance (15 by 15 cells, time horizon 16, and three searchers), the resulting algorithm 
essentially provides an approximate solution that is guaranteed to be within 5% of an 
optimal solution in less than about 20 minutes. Instances with more searchers (such 
as 10, 15 and 30 searchers) are solved even faster. In addition, we prove that under 
certain assumptions the nonlinear convex multiple homogeneous searcher problem is 
equivalent to a large-scale linear mixed-integer program. 



I. INTRODUCTION 

A. MOTIVATION 

The need to search for moving objects arises in many civilian and military 
applications. Rescue teams search for lost persons and shipwrecks. Autonomous 
robots could hud victims and survivors after natural disasters. In military search 
operations, patrol aircraft and unmanned aerial vehicles (UAVs) look for high-valued 
targets such as missile launchers and terrorist vehicles. In all these applications, 
the choice of paths for the searchers strongly influences the probability of hnding 
the targets within a specihc time limit. Unfortunately, the problem of selecting the 
“best” path is fundamentally hard due to the nonlinearity induced by the probability 
of detection. For example, looking twice by a searcher generally does not double the 
detection probability. 

In military search missions with UAVs, the problem of finding good paths 
for searchers (UAVs) is complicated even further. Rather than manned aircraft, 
UAVs are preferred to operate in dangerous environments such as hostile regions and 
battlehelds. In these applications, the searchers (UAVs) often need to balance search 
effectiveness with threats to themselves posed by surface-to-air missiles and small- 
arms hre. Moreover, the search problem becomes signihcantly more difficult as the 
number of searchers grows. With the technological advances in cooperative control of 
multiple UAVs, such multi-searcher problems are appearing with increasing frequency 
in applications. The need for optimal path planning in these complicated situations 
motivates this dissertation. 

In this dissertation, we consider several search problems including situations 
with searchers subject to threats as well as multiple searchers. We particularly focus 
on search by aerial assets that are subject to threats from the ground. We formulate 
search problems for single and multiple searchers and develop solution methodologies 
for efficiently hnding an optimal or near-optimal path. 
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B. SCOPE OF DISSERTATION 

We consider discrete time-and-space path-optimization problems where a non¬ 
evading target moves in a two-dimensional area of interest (AOI). The AOI is dis¬ 
cretized into a hnite set of cells and the target moves between the cells in discrete 
time. Stone [47] and Brown [7] introduced the basic forms of this search problem. 
A single searcher or multiple searchers move, in discrete time, through a discretized 
two-dimensional space (on the ground in the AOI or at a constant altitude over the 
AOI) or a discretized three-dimensional airspace over the AOI. The searchers’ space 
is represented by a directed graph, where vertices represent waypoints (or search 
sectors) from which a searcher can search a particular cell in the AOI and where di¬ 
rected edges represent transition between two waypoints. Thus a path of the searcher 
is represented as a sequence of waypoints (vertices). 

We note that our aim is not to hnd an optimal flyable trajectory that fully 
accounts for turn radius, climb/dive angle, and max/min speed for the searcher. We 
aim for “high-level” control of the searchers where waypoints may represent points 
in space that are many minutes apart in flying time and may represent areas to be 
searched for some period of time. We consider imperfect searchers that may search a 
cell that contains a target without hnding the target. However, in this dissertation, 
we assume that the searchers will not report a target in a cell where there is no target. 

The goal for the searchers is to maximize the probability of detecting the target 
within a hnite mission time or, equivalently, to minimize the probability that the 
target is not detected during the duration of the mission. The searchers are subject 
to constraints on path continuity, mission time, and possibly other “resources” such 
as risk exposure and fuel consumption. 

This dissertation refers to the path optimization problem for a single searcher 
with path- and time-constraints as the single searcher problem (SSP). Other authors 
refer to this problem as the optimal searcher path problem [41, 49, 32]. In this 
dissertation, we focus on the following variants of SSP. 
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1. SSP with additional constraints on, e.g., risk exposure to threats and fuel 
consumption during the missions. We refer to this path-optimization problem 
as the resource-constrained single searcher problem (RSSP). 

2 . SSP for multiple searchers where the probability of detecting the target by at 
least one searcher is maximized. We term this problem the multiple-searcher 
problem (MSP). 

3. MSP for homogeneous searchers where the fact that searchers are identical is 
utilized. We refer to this multiple-homogeneous-searcher problem as MHSP. 

This dissertation always assumes that the searchers are subject to constraints on 
path continuity and mission time. We note that the multiple-searcher problem with 
resource-constraints (risk exposure and fuel consumption, etc.) is outside of the scope 
of this dissertation. 

C. LITERATURE SURVEY 

The topic of searching for a moving target in discrete time and space where 
a searcher’s path is constrained has received much attention. Benkoski et al. [41] 
provide a comprehensive survey of this problem (until 1990). The path- and time- 
constrained search problem with a stationary target is NP-complete [49]. Thus the 
path- and time-constrained, moving-target, search problem is at least as difficult. 

One scheme to hud an optimal path for a searcher is the branch-and-bound 
procedure. A branch-and-bound procedure was hrst studied by Stewart [46]. Since 
Stewart’s “bounds” are not valid, an optimal branch of the enumeration tree may be 
fathomed. More recent studies [18, 34, 52, 32] focus on the development of specialized 
branch-and-bound algorithms for hnding an optimal path for the single searcher. In 
these algorithms, a path is a sequence of vertices that the searcher will visit and 
branching corresponds to extending a subpath by one more vertex. Bounds on the 
optimal value of the problem are obtained by replacing the probability of detection 
with, effectively, the expected number of detections [34, 52, 32]. Bounds are also 
obtained by assuming that the searcher can divide its effort among multiple cells each 
time period [18]. The efficiency of various branch-and-bound methods is investigated 


3 



in [52], with emphasis on the tradeoff between the accuracy of the bound employed 
and the time required to compute it. 

An alternative approach to obtain an optimal path is based on dynamic pro¬ 
gramming and partially observable Markov decision processes [17]. The approach 
can solve small problems quickly, but requires a large amount of computer storage as 
the problem size increases. Later the same author reports that a branch-and-bound 
algorithm [18] performs more efficiently and can solve larger problems than the exist¬ 
ing dynamic programming procedures. Thus, we mainly focus on branch-and-bound 
procedures in this dissertation. 

1. Resource-constrained Single Searcher Problem 

In the studies listed above, it is assumed that the searcher moves on the 
ground in the AOI or searches the AOI from a constant altitude. However, in mil¬ 
itary search, surveillance, and reconnaissance operations over land, factors such as 
risk and fuel consumption become independent elements of concern for planners [37], 
and the searcher is required to change its flight altitude during a mission. Risk to 
the searcher arises from exposure to enemy threats on the ground such as small-arms 
hre, anti-aircraft artillery, and surface-to-air missiles. Risk of pilot error (e.g., for a 
low-flying helicopter [29]) and mechanical failure (e.g., for small and unreliable UAVs 
[31]) may also be significant. Fuel consumption tends to be proportional to the search 
duration. However, for small UAVs, it is not always proportional to the mission time, 
due to variation in the searcher’s speed and/or altitude, as well as varying weather 
conditions. During search over land, terrain features require altitude changes. The 
altitude is also varied to balance the searcher’s risk with image quality. We note 
that image quality is of particular concern for small UAVs operating with low- to 
moderate-quality sensors (For information about UAV surveillance operations and 
sensors, see [37, 42]). In civilian search and surveillance operations over land, we may 
face many of the same factors with the exception of ground fire. 
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There are only a few studies dealing with constraints on risk and fuel con¬ 
sumption for the searchers. Thomas and Eagle [48] consider a search problem with 
a random time horizon; i.e., with a perishable target with a random lifetime that 
follows a particular distribution (a geometric distribution in that study). This situa¬ 
tion is similar to having the searcher being subject to threats and having the search 
terminate after searcher “failure.” We will consider a similar situation, but will limit 
the total risk a searcher can be exposed to during the mission. Secrest [43] considers 
optimal path planning with multiple-objectives (simultaneously accounting for flight 
time, risk and fuel consumption, etc.) in surveillance missions while visiting all as¬ 
signed locations. The resulting model does not consider the probability of detecting 
a target. 

Models for military aircraft routing frequently consider minimum risk along a 
path subject to fuel consumption and mission duration constraints, see, e.g., [35, 55, 
10]. These models have similar constraints to the ones in RSSP, but deal with linear 
and time-invariant objective functions. In contrast, RSSP has a nonlinear, time- 
variant objective function representing the probability of detection (see Chapter II). 
Only recently have researchers extended the models for aircraft routing to situations 
with nonlinear objective functions, and then only for special cases of path-dependent 
risks [29]. 

2. Multiple-Searcher Problem 

The branch-and-bound procedure is also applicable in solving the multiple- 
searcher problems (MSP/MHSP). Dell et ah [13] present an exact procedure (branch- 
and-bound algorithm) and six heuristics (local search, expected detection heuristics, 
genetic algorithms, and moving time-horizon heuristic) to solve the MSP/MHSP. 
However, the branch-and-bound algorithms need prohibitive runtime to guarantee an 
optimal path for multiple-searchers, because the number of possible paths grows expo¬ 
nentially with the number of searchers. This motivates the development of heuristic 
algorithms. We find other heuristic algorithms based on myopic and receding-horizon 
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approaches in [23] and sequential optimization of each searcher in [26]. For a specihc 
scenario, Riehl et ah [28] present a receding-horizon approach that jointly optimizes 
paths and sensor orientations for multiple searchers. For a stationary-target case, 
Song et ah [45] use a forward induction approach to hnd an optimal search strategy. 

Control and robotic communities have analyzed the search problems for mul¬ 
tiple UAVs and sensors as well. Ryan et al. [2] review the current issues and develop¬ 
ments in the held of cooperative control (real-time navigation, collision avoidance, and 
hight formation, etc.). Decentralized search techniques [20, 54] and Bayesian heuristic 
approaches [53, 25] are also studied in the context of real-time search operations. 

D. CONTRIBUTIONS 

This study considers the path optimization problems, including “searcher sub¬ 
ject to threat” and cases with multiple searchers. We extend the classical single 
searcher problem (SSP) to realistic situations such as: (1) the searcher needs to 
change its hight altitude during a mission to balance sensor quality and risk from the 
ground treats; (2) the searcher is subject to constraints on multiple resources such as 
fuel consumption and risk exposure along its path; and (3) the search ehect depends 
on the current and previous locations of the searcher. The latter situation may occur 
if moving from some search sector to a new sector results in a lower search ehect due 
to inehective search during transit. 

We present new bounding techniques and network reduction procedures in a 
branch-and-bound algorithm to efficiently solve the SSP. Based on this algorithm, we 
develop a specialized branch-and-bound procedure to solve the resource-constrained 
single searcher problem (RSSP) by utilizing several novel network reduction proce¬ 
dures as well as new bounding technique taking account of the multiple resource 
constraints. We also develop several new algorithms to solve the multiple searchers 
problems (MSP/MHSP). 
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E. ORGANIZATION 

The reminder of this dissertation is outlined as follows. In the next chapter, 
we formulate RSSP and develop a specialized branch-and-bound algorithm for RSSP. 
The resulting algorithm utilizes a new bounding technique, applies several network 
reduction procedures, and exhibits promising behavior in computational tests on in¬ 
stances of SSP and RSSP. In Chapter III, we formulate MSP by extending SSP and 
present three new algorithms (an exact procedure and two heuristics) to solve MSP 
and examine their computational efficiency. In Chapter IV, we focus on MHSP. We 
use the convexity of the nonlinear objective function (the non-detection probability) 
and construct an exact algorithm using cutting planes (outer approximations). We 
prove that under certain assumptions the problem is equivalent to a large-scale lin¬ 
ear mixed-integer program. We also introduce several new cuts for the MHSP and 
demonstrate their effect in computational tests. Chapter V provides conclusions of 
our study as well as suggestions for future work. 
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II. RESOURCE-CONSTRAINED SINGLE 
SEARCHER PROBLEM 

This chapter formulates the resource-constrained single searcher problem (RSSP) 
by generalizing existing models of the single searcher problem (SSP) to the case 
with multiple resource constraints. The RSSP also considers an altitude-dependent 
searcher with its search effect depending on not only the searcher’s current loca¬ 
tion but also its previous location. We combine the branch-and-bound methodology 
for solving SSP with the line of research on the resource-constrained, shortest-path 
problem [10]. Specihcally, we merge the algorithms in [32] and [10], and develop a spe¬ 
cialized branch-and-bound algorithm for RSSP. The resulting algorithm utilizes new 
bounding techniques, applies several new network reduction procedures, and exhibits 
promising behavior in computational tests on instances of SSP and RSSP. 

A. PROBLEM DESCRIPTION AND FORMULATION 

The area of interest (AOI) is discretized into a hnite set of cells C = {1,..., C*} 
(see Figure 1) and the time horizon is discretized into a hnite set of time periods 
T = {1,2, ...,T}. A target occupies one cell each time period and moves among 
cells according to a Markov process with known transition matrix P = (Pc,c')) where 
Pc,c' is the probability that the target moves from cell c to cell c' dnring one time 
period. This Markovian target motion is hrst modeled in Stone [47]. Let p{-,t) = 
[p{l,t),p{2,t),... ,p{C,t)], where p{c,t) is the probability that the target is in cell 
c G C at the beginning of time period t E T and the target has not been detected 
before t. We refer to p{-,t) as the undetected target distribution. We note that p{-,t) 
may not be a probability mass fnnction as the snm of the elements could be less than 
1. The initial target distribntion p{-, 1) is known. 

A single searcher moves through a designated airspace over the area of interest 
with the goal of hnding the moving target on the gronnd. The airspace over each 
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Figure 1. A discretized area of interest composing of C* = 4 cells. 

cell c G C is vertically discretized into a set of altitudes (or boxes) Td = {1, 2,..., H} 
(see Figure 2). For any c G C and h eH, we refer to the cell-altitude pair (c, h) as a 
waypoint where the searcher can loiter and carry out search of cell c. We assume that, 
from waypoint (c, h) , the searcher can sweep cell c only. We model the designated 
airspace by a directed network (V,£^), with set of vertices V and set of directed 
edges S, in which vertices v = {c,h) E V represent waypoints and directed edges 
e = (u, v') E S represent transition between waypoints v, v' E V. The waypoints 
v,v' eV are adjacent to each other if v and v' touch at any portions. The searcher 
can only transit between two waypoints that are adjacent to each other. Let J^{v) eV 
be the set of vertices that are adjacent to u G V. We refer to J^{v) as the forward 
star of vertex v. We adopt the convention that v E J^{v) for all v eV. Then, the set 
of edges £ = {(u,u') G V x V| F G 

During each time period t E T, the searcher is at a particular vertex (way- 
point). We assume there is no transit time between waypoints. Hence, {v,v') E £ 
simply represents search at waypoint v followed by search at waypoint v' in the next 
time period. The situation with nonzero transit time between waypoints can be mod¬ 
eled, at least approximately, by introducing artificial vertices. (We refer to [32] for 
a comprehensive study of nonzero transit times.) We note that the edge {v,v) E £ 
represents searching at waypoint v for two consecutive time periods. 
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Figure 2. A discretized airspace over the area of interest (Figure 1) with H = 2 
altitudes. 

Let 0 : V —>■ C be the function that specifies the cell over which a vertex is 
located, i.e., cell 0(u) is searched from vertex v. We denote the searcher’s vertex 
(waypoint) prior to time period 1 by uq ^ V. This starting vertex could be a desig¬ 
nated entry waypoint into the area of interest. We also dehne V C V to be a set of 
possible destination vertices for the searcher (e.g., V = {uq} if the searcher is required 
to return to the starting vertex or V = V if the search can end anywhere). 

For any k E T and Vi E V, I = 0,1,2,k, such that {vi-i,vi) E £ for all 
I = 1,2,...,k, let the sequence denote a directed v^-Vk subpath. If Vk E V, 

then the directed VQ-Vk subpath is a directed Vo-Vk path. When the meaning is clear 
from the context, we refer to a directed vo-Vk (sub)path as a (sub)path. In this 
notation, the searcher flies from vq to some Vk E V along a directed Vo-Vk path. The 
searcher occupies only one vertex v E V each time period, and stays at the same 
vertex or moves to another vertex in lF{v) for the next time period. 

We adopt the following target-detection model. If the target is in cell c E C 
during time period t E T and the searcher is at the same time at waypoint v' E V 
above cell c, i.e., (/>(u') = c, then detection occurs with a probability. We term this 
probability as glimpse detection probability and refer to g{v,v',t), where u G V is the 
searcher’s waypoint during time period f —1. Hence, the glimpse detection probability 
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during time period t depends on the previous and current waypoints for the searcher. 
This is a generalization of earlier models where the glimpse detection probability 
depends only on v' and t [52, 32] and is one of our contributions. Our model accounts 
for the fact that moving from some waypoint to a new waypoint may result in a lower 
glimpse detection probability than if the searcher already was loitering at the latter 
waypoint, i.e., g{v',v',t) > g{v,v',t) for v' ^ v. This effect occurs if refocusing the 
sensor and becoming familiar with a new cell have a signihcant detrimental effect on 
the glimpse detection probability. In general, change of waypoint, especially change 
of altitude and frequent, irregular change of direction, may distract from the search. 
This generalization also allows us to account indirectly for small transit times (much 
less than the length of a time period) between waypoints without adopting a hne time 
discretization with resulting high computational cost. In this notation, the probability 
of detection at waypoint v' during time period t, given search at waypoint v during 
time period t — 1, and no prior detections becomes p{(f){v'),t)g{v,v',t). 

The glimpse detection probability may also depend on the searcher’s speed. 
To account for this situation, edges can be duplicated as in [10] to represent search 
at different speeds in cases where speed influences the searcher’s effectiveness signif¬ 
icantly. For the sake of simplicity, however, we do not introduce notation to handle 
that situation. 

Since p{-,t) is the undetected target distribution at the beginning of time 
period t, it depends on searches up to time period t — 1. Specihcally, if p{-,t) = 
[p{l,t),... ,p{c',t),... ,p{C,t)] and cell c' is searched from waypoint n' (i.e., 0(n') = c') 
during time period t, the undetected target distribution at the beginning of the next 
time period t -|- 1 is 

p{-,t + l) = [p{l,t),..,p{c-l,t),p{c,t){l-g{v,v',t)),p{c+l,t),..,p{C,t)]r, (II.l) 
where v is the searcher’s vertex during time period t — 1. Thus, the probability of 
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detection along a directed Vo-Vk path V = denoted q(V), is defined as 

k 

Qi'P) = ^P{<P{vt),t)g{vt-i,vt,t). (II.2) 

t=i 

Let X = {1,2,...,/} be a set of resources and fi{v,v\t) be the amount of 
resource i E I consumed by the searcher at vertex v' during time period t, given search 
at vertex v during time period t — 1. Resources may represent physical commodities 
such as fuel and ordnance that are depleted during the search as well as abstract 
factors such as notions of risk exposure during the search. In contrast to search by 
manned aircraft, where significant risks are usually avoided, planners accept higher 
risks for UAV search missions and would like to balance risk with other factors during 
the planning process. We discuss the calculation of risk in Section D. The total 
“consumption” of resource i E I along the directed path V = {uzjjLo 

t). (II.3) 

t=i 

The searcher cannot consume more than fj of resource i E X along a path. Hence, 
the resource-constrained single searcher problem (RSSP) is the problem to find a 
directed path V = (u/ULq, with Vk eV, k E T, that maximizes qiV) subject to 
the constraints 

TiiV) <fi,i eX. (II.4) 

We contribute that the single searcher problem (SSP) is generalized to the case 
with multiple resource constraints. We refer to constraints (II.4) as side constraints. 
In Section D, we examine a case study with two time-invariant resources: risk and fuel. 
The next section deals with the “unconstrained” problem with no side constraints 
(referred as SSP), while Sections C and D address the full RSSP. 

B. BRANCH-AND-BOUND ALGORITHM FOR SSP 

In this section, we consider the single search problem (SSP) that maximizes 
(II.2) without the side constraints (II.4). Several branch-and-bound algorithms for 
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the solution of SSP have been developed [46, 50, 18, 34, 13, 52, 32], These algorithms 
implicitly enumerate all searcher paths that cannot be proven, by means of a bound, to 
be nonimproving. The next subsection presents the algorithm in [32], which appears 
to be the fastest in the literature for SSP under a broad range of conditions; see [32] 
for an empirical study. The subsequent snbsection presents fonr modihcations of that 
algorithm which generalize and speed up that algorithm. 

This section ignores side constraints and, without loss of generality, assnmes 
that an optimal path consists of T + 1 vertices. To simplify the notation, we also 
assume in this section that there is no end-point restriction, i.e., V = V. We hnd 
identical assumptions in [52, 32]. Initially, we assnme that the glimpse detection 
probability g{v,v',t) is independent of v and write g{v',t), but relax that assumption 
later in this section. 

1. Existing Algorithm 

For completeness and ease of reference later, we ontline the algorithm in Lau 
et ah [32]. Given a subpath {viYiZq, t E T, let /C(t) be the set of triplets of the form 
q{vt,t)) representing extensions of yet to be explored. The first element 

Vf refers to the next vertex to visit, the second element t is the time period^ to visit 
the vertex Vt, and the third element q{vt,t) is an npper bound on the probability 
of detection along any path that starts with the subpath {vi}Yo- The upper bound 
qivtit) consists of three parts. Let dt{vt,t) be an upper bound on the probability of 
detection during time periods t-|-l,t-|-2,...,T, given that the searcher starts at Vt at 
time t, and no detection occnrs along the snbpath {vi}Yo- The two other parts are 
the probability of detection on the subpath and the probability of detection 

dnring t. Hence, 

q{vt,t) = q{{vi}\zl) + p{(t){vt),t)g{vt,t) + dt{vt,t). (II.5) 

^This information is currently redundant but the notation is convenient in later generalizations. 
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We also let q denote the largest detection probability found so far among all the 
examined paths. In this notation, the algorithm in [32] takes the following form. 

Algorithm 1 (Solves SSP). 

Input: A time-expanded network {Af,A) (dehned later as Figure 3) with 
nodes n E Af and arcs (n,n') G A where n = {v,t — 1), n' = 

Arc lengths g{v', t) in (A/”, A), the initial target distribution p(-, 1), Markov 
transition matrix F, and upper bound q. 

Output: An optimal search path V* = {vi}f^Q and its value q*. 

Step 0. Set t = 0, IC{t) = {(uq, 0, cxd)}, and q = 0. 

Step 1. If )C(t) is empty, replace t hj t — 1. Else, go to Step 3. 

Step 2. If t = 0, stop: the last saved path is optimal and q is its probability 
of detection. Else, go to Step 1. 

Step 3. Remove from IC(t) the triplet {vt,t,q{vt,t)) with the largest bound 
Q{vut). 


Step 4. If q{vt,t) < q, go to Step 1. (Current subpath is fathomed.) 

Step 5. If t < T, then for each vertex v G Alvt), calculate a bound dt+i(v,t + 

1) as well as q(v, t -|- 1) using equation (II.5), and add (v, t + 1, q{v, t -|- 1)) 
to Kit + 1). Replace thy t + 1 and go to Step 3. Else, let q = q{vt, t) and 
save the incumbent path {vi}JLq, and go to Step 1. 

Clearly, a tight bound dt{vt,t) will reduce the number of branching attempts 
in Step 4 of Algorithm 1. As examined in [51, 52, 32], there is a fundamental trade-off 
between the effort needed to compute a bound and its tightness. From these studies, 
it appears that the bounding technique in [32] compares favorably in most situations. 
We describe that bounding technique in the rest of this subsection. 

Consider a subpath {viYi^q, t E T, and let Pg{-,t) be the undetected target 
distribution after search along {vi}Yo, i-®-) 

Pg{-,t)= (II.6) 

[p(l,t),..,p(0(nt) - l,t),p{(f){vt),t){l - g{vt,t)),p{(l){vt) + l,t),..,p{C,t)]. 
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We use subscript g to indicate that Pg{-,t) is obtained from p{-,t) by applying the 
glimpse detection probability corresponding to the last vertex in the “current” sub¬ 
path {viYi^Q. For any integer s > t, s E T, we also dehne 

pr{;s;t)=Pg{;t)r-\ (II.7) 

As seen, pr{c,s-,t) is the probability that the target is in cell c at time period s > t 
and there was no detection during search along the subpath {vi}Yo- Hence, target 
distribution pr{-,s]t) at time period s > t ignores the effect of search after time 
period t. If the subpath is {uq}, i.e., t = 0, we dehne for notational convenience 

pr(-,5;0)=p(-,l)r-\ (II.8) 

for any s > 0, s G T. Moreover, we dehne pr(c, t;t) = 0 for all c G C and t = 0,1,..., T. 

Now, we construct a time-expanded graph from the network (V,£^) as follows 
(see Figure 3). Each vertex u G V is duplicated T times to dehne the nodes {u,<s), 
s E T. Let Af be the set of all such nodes as well as the nodes uq = {vq, 0) and 
n = {v,T+l) representing the searcher’s prior position and hnal position, respectively. 
Here, v is an artihcial terminal vertex. Two nodes n = {v,s — 1) and n' = {v',s), 
v,v' eV and s = 2, 3,..., T, are connected with an arc (n, n') if and only if (u, v') E S. 
Moreover, the node uq = (uo,0) is connected with an arc to a node n' = {v',1), 
v' E V, if and only if (vo,v') E and every node n = (u,T), u G V is connected 
with an arc to fi. Let A be the set of all arcs. For any integer k <T + 1 and nodes 
ni = {vi, 1) E Af, / = 0,1,..., k, such that (n^-i, ni) E A for all I = 1,2, ...,k, we let the 
sequence denote a subpath in the time-expanded graph {Af,A). 

For some t E {0,1, ...,T — 1}, suppose that a subpath {vi}Yo original 

graph {y,£) is given. Then, we endow each arc {n,n') = {{v,s),{v',s + 1)) G A, 
s = t,t + l,...,T — 1, in the time-expanded graph {Af, A) with a “reward” 

Cn,n' = [PT{(t>{v'), s + l-,t) - s; t)g{v, s)F(u, v')]g{v', s -f 1), (11.9) 

where F(u,u') is the (j){v)-(j){v') element of the Markov transition matrix F. We set 
Cn,h = 0 for all (n, n) E A. We observe that, multiplied out, the hrst term pv{(t>{v'), s-f 
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Figure 3. A time-expanded graph from the network in Figure 1 and one altitude. The 
searcher’s initial position is Uq = (uq, 0) = ((1,1), 0) and hnal position is h = (h, T-l-1). 

l]t)g{v',s -|- 1) in (II.9) is effectively the expected number of detections during time 
period s -|- 1, which gives rise to the so-called mean bound [34], and the second term 
in (II.9) improves the bound by accounting for the effect of search during time period 
s [32]. We refer to {N,A) with arc rewards given by (II.9) as the time-expanded 
network. 

In Lau et al. [32], it is shown that given the subpath the optimal value 

of the longest-path problem in the time-expanded network from node (ut,t) to node 
('0,T -|- 1), using the rewards in (II.9) as “arc length,” provides an upper bound for 
Algorithm 1. Specihcally, this optimal value is an upper bound on the probability of 
detection during time periods f-|-l,t-|-2,...,T, given that the searcher starts at Vt at 
time t, and no prior detections occurred along the subpath We denote this 

bound by dt{vt, t) and refer to it as the dynamic bound, as it needs to be recomputed 
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every time the current subpath is extended. 


Since the time-expanded network is acyclic, the longest-path problem can be 
solved in polynomial calculation time with a standard shortest-path algorithm ([1], 
pages 77-79). We observe that to compute dt{vt,t) given the subpath it is 

only necessary to generate the part of the graph {J\f,A), and the corresponding arc 
rewards, “after” time t and within reach from node {vt,t), since the longest path 
starts at node {vt,t)- Hence, in Step 5 of Algorithm 1, it suffices to generate the 
time-expanded network “after” time t. 

2. Algorithmic Modifications for SSP 

Algorithm 1 requires bound calculation at each branching which is a principal 
bottleneck for that algorithm. Thus we present new bounding techniques to enhance 
Algorithm 1. This subsection proposes and examines four modifications of Algorithm 
1. The hrst modification extends the bound in [32] to the case in which glimpse 
detection probability depends on the previous vertex, i.e., g{v,v',t). The following 
three modihcations are aimed to simplify and speed up Algorithm 1. The second 
modihcation weakens the bound but greatly simplihes its bound calculation. The 
third modihcation improves the weaker bound at little computational expense. The 
fourth modihcation takes advantage of a special, but frequently occurring, initial tar¬ 
get distribution. 


a. Bound for Edge-Dependent Glimpse Detection Proba¬ 
bility 

Previous studies (see, e.g.,[52, 32]) assume that the glimpse detection 
probability depends only on the searcher’s current waypoint (vertex) and on time. 
As we argue in Section A, this is somewhat restrictive. Fortunately, we can easily 
extend the previous studies to the case where the glimpse detection probability also 
depends on the searcher’s previous waypoint. The only modihcation that is required 
is to redehne the arc reward Cny in (II.9). However, a straightforward replacement of 
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g{v, s) by giy”, v, s) in (II.9), where v” is the vertex prior to v, would ruin the longest- 
path structure of the bound-calculation problem; Cny would no longer depend only on 
the head and tail of the arc (n, n'). Hence, it is necessary to use the smallest glimpse 
detection probability min„//g 7 ^(„) g(v", v, s) to eliminate the dependence on the vertex 
prior to v, where 7^(w) C V is the reverse star of v, i.e., 7^(v) = {v" G V | (u", v) E £}. 
Consequently, we now endow each arc (n, n') = ((w, s), {v\ s-|-l)) G A with the reward 


Cnri/ — 


Pr(0(v),s + l;t) -pr{(l){v),s;t) min g{v”,v,s) T{v,v) 

\v"e:TZ{v) / 


g{v,v',s + l). 

( 11 , 10 ) 

With this reward, the bound calculation remains a longest-path problem in an acyclic 
graph and it can be shown using the same arguments as in [32] which prove that the 
dynamic bound is valid. 


b. Static Bound 

Algorithm 1 requires one longest-path calculation in the time-expanded 
network for each vertex in the forward star of the current vertex to compute the re¬ 
quired bounds dtiytit) (see Step 5). The approach of Algorithm 1 with dynamic 
bound follows the traditional approach of branch-and-bound algorithms where the 
bound is reoptimized before each branching. In the present case, the reoptimization 
corresponds to the longest-path calculation and requires computing the arc rewards 
Cn,nG see (II. 10). These calculation can be time consuming, as they typically involve 
a moderately large Markov transition matrix T and associated matrix multiplication. 
Our numerical example considers that T is of size 225 by 225. We propose to use a 
static bound instead of the dynamic bound proposed in [32] and described in Subsec¬ 
tion B.l. As shown below, all the necessary static bounds are computed prior to any 
branching and are not recomputed later. 

The dynamic bound dt{vt,t) (see Subsection B.l) uses information 
about search along a current subpath {viYi^q. However, the bound remains valid 
if the effect of search along the current subpath is ignored. This follows from the 
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same arguments as in the proof of the validity of dt{vt,t), see [32], We denote the 
new bound do{vt,t), where the subscript 0 indicates that the trivial subpath {uq} 
is used in (II.10) with t = 0 instead of the subpath {viYi^q, which dt{vt,t) utilizes. 
Hence, the only difference between dt{vt, t) and d^ivt) t) is that pr{-, ■; t) is replaced by 
Pr(',s0) in (II.10). However, this difference makes the bound do{vt,t) independent 
of the current subpath used to reach the vertex u*. Hence, do{v,t) can be computed 
in advance for all nodes {v,t) G Af, and dynamical computation of bounds is not 
required. Consequently, the arc rewards (11.10) and bounds are computed only once. 
We refer to do{v,t) as the static bound. We observe that it is not necessary to carry 
out a longest-path calculation from each node {v, t) G A/” to (h, T + 1) to obtain 
do{v,t). It is more efficient to carry out the longest-path calculations backward from 
node (h, T 4-1) to all nodes. This calculation simply amounts to applying a shortest- 
path algorithm once to the time-expanded network with arc lengths equal to the 
negative rewards. 

In Step 5 of Algorithm 1, we now simply use do{vt, t) instead of dt{vt, t). 
Thus, the modified algorithm does not require any longest-path calculation in Step 5. 
All bound calculations are done prior to Step 0. Clearly, the modified approach results 
in a weaker bound and the need for more branching attempts. However, the additional 
branching attempts may be compensated by shorter per-iteration computing times. 

In order to examine the effect of the static bound do{vt,t), we examine 
the same numerical example as in [32]: An area of interest consists of 11 by 11 cells. 
The searcher operates only at one altitude and its moves are restricted to vertically 
and horizontally adjacent cells, excluding diagonal moves. The target remains in the 
current cell with a probability p or moves to one of the vertically or horizontally 
adjacent cells with probability 1 — p. The different moves are equally likely. The 
searcher departs cell 1 (uq = 1), where the cells are numbered from left to right and 
from top to bottom. Hence, cell 1 is the upper-left-corner cell. The target starts at 
the center cell, i.e., at cell 61. The time horizon T = 17. 
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We have implemented Algorithm 1 with static bound using Microsoft 
Visual C++ 6.0 on a desktop computer with a 3.4 GHz Intel Pentium IV processor, 
1.0 gigabytes of RAM, and the Microsoft Windows XP operating system. Table 1 
shows, for a range of constant glimpse detection probabilities g{v,v',t) and “stay 
probabilities” p, the run times (in seconds) and numbers of bounding attempts of 
Algorithm 1, as well as those reported in [32]. We especially consider the case of high 
glimpse detection probability {g{v,v',t) = 0.99), in which both static and dynamic 
bounds tend to be relatively weak. Column 3 of Table 1 presents the resulting run 
times for Algorithm 1 with static bound. The corresponding run times reported from 
[32] are found in column 5. (The case g{v,v',t) = 0.99 is not considered in [32].) 
Those reported numbers are achieved on a 2.6 GHz Opteron 152 processor computer 
using Matlab. Hence, a direct comparison between the run times in columns 3 and 
5 is difficult. However, we hnd reports of run times for other problem instances 
using a C++ implementation in [32]. A brief comment in [32] based on a single 
problem instance indicates that the Matlab implementation is 15.6 times slower than 
the C++ implementation. For the sake of comparison, we scale down the run times 
in column 5 of Table 1 with 1/15.6 to approximately account for the slower Matlab 
implementation. We also scale down the run times of column 5 by a factor of 2.6/3.4 
to account for the slower computer used in [32]. The resulting scaled down run times 
are presented in column 6 of Table 1. We observe that the scaled down run times for 
Algorithm 1 with dynamic bound are somewhat faster than the corresponding run 
times with static bound. However, the static bound, which is weaker than the dynamic 
bound, remains fairly competitive, especially for more difficult problem instances. 

We compare the strengths of the static and dynamic bounds by counting 
the number of branching attempts required in Algorithm 1. Columns 4 and 7 of 
Table 1 give the numbers of branching attempts for the static and dynamic bounds, 
respectively. The number of branching attempts for the static bound is, on average, 
37 times larger than in the case of the dynamic bound. We observe that the greater 
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Static Bound 

Dynamic Bound [32] 

Dynamic Bound 



Time 

Branching 

Time 

Scaled 

Branching 

Time 

Branching 

g{v,v',t) 

P 

(sec.) 

attempts 

(sec.) 

(sec.) 

attempts 

(sec.) 

attempts 


0.3 

3.59 

2,301,182 

14.56 

0.71 

58,314 

2.09 

58,314 

0.3 

0.6 

2.58 

1,635,517 

14.53 

0.71 

52,369 

2.11 

52,369 


0.9 

11.47 

7,338,492 

157.30 

7.71 

380,889 

20.75 

380,889 


0.3 

5.38 

3,424,282 

19.07 

0.94 

49,774 

2.59 

49,774 

0.6 

0.6 

2.58 

1,620,402 

23.76 

1.16 

47,454 

3.03 

47,454 


0.9 

59.26 

38,186,809 

730.96 

35.84 

2,185,066 

103.08 

2,185,066 


0.3 

5.78 

3,675,197 

21.55 

1.06 

45,019 

2.78 

45,019 

0.9 

0.6 

2.89 

1,831,875 

30.58 

1.50 

59,527 

3.91 

59,527 


0.9 

176.26 

113,646,357 

2902.27 

142.30 

11,299,259 

431.38 

11,299,259 


0.3 

6.09 

3,865,002 

N/A 

N/A 

N/A 

2.91 

46,339 

0.99 

0.6 

3.00 

1,896,960 

N/A 

N/A 

N/A 

3.91 

58,752 


0.9 

192.06 

123,822,672 

N/A 

N/A 

N/A 

558.63 

16,685,969 


Table 1. Run times and number of branching attempts (counted in Step 4) for 
Algorithm 1 with static and dynamic bounds on 11 by 11 cell search problem with 
time horizon T = 17. Columns labeled “Dynamic Bound [32]” correspond to original 
and speed-adjusted results from [32], 


number of bounding attempts is partially compensated for by avoiding dynamical 
reoptimization of the bound. 

For a direct comparison between the static and dynamic bounds, we 
implement Algorithm 1 with dynamic bound in Microsoft Visual C-I--I- 6.0. We made 
a signihcant effort to ensure that the implementation is efficient, including efficient 
handling of sparse matrices. Column 8 and 9 of Table 1 report the run times and the 
number of branching attempts for our implementation of Algorithm 1 with dynamic 
bound, respectively. We observe that our implementation results in identical numbers 
of branching attempts compared to the implementation in [32]. When comparing 
columns 6 and 8, we see that our implementation of the dynamic bound results in 
somewhat longer run times than the scaled times from [32]. However, the longer times 
in column 8 compared to column 6 can partially be due to an excessively aggressive 
scaling of run times going from column 5 to column 6. 

While implementing the dynamic bound, we noted a signihcant chal¬ 
lenge associated with efficient matrix multiplication and data handling. Since the 
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dynamic bound dt{v,t) is reoptimized a large number of times, it is paramount to 
carry out the related calculations efficiently. In comparison, it is rather trivial to 
implement a static bound. It is not critical to carry out the longest-path calculations 
highly efficiently in the static case as they are only done once. 

We also observe that both static and dynamic bounds tend to be weaker 
when the target is nearly stationary (e.g., p = 0.9) and glimpse detection probability 
is large (e.g., g{v,v',t) = 0.9 or 0.99). For these problem instances, the implementa¬ 
tion with dynamic bound is more sensitive to the change in data. In fact, comparing 
the instance {g{v,v',t) = 0.9 and p = 0.9) to the instance {g{v,v',t) = 0.99 and 
p = 0.9), we see that the run time and the number of branching attempts in the 
case of the dynamic bound become 1.29 and 1.48 times larger, respectively. On the 
other hand, the static bound is less sensitive, and its numbers become only 1.09 times 
larger. These numbers indicate that it becomes even less worthwhile to invest time 
in dynamic bound calculations when the bounds are relatively weak anyway. 

c. Directional Static Bound 

As seen from Table 1, the static bound is substantially weaker than 
the dynamic bound. We derive a stronger static bound motivated by the classical 
approach to handling turn-radius constraints in vehicle-routing problems [8]. 

In the longest-path calculations for the static bound, the reward of arc 
{{v,s), {v',s + 1)) is, effectively, the probability of detection at vertex n' during time 
period s -|- 1 and no detection at vertex v during time period s. Of course, this 
overestimates the probability of detection at vertex v' during time period s -|- 1 and 
no prior detections, as detection could occur prior to time period s. We strengthen 
the static bound if we redehne the arc reward to be the probability of detection at 
vertex v' during time period s -|- 1 and no detection at vertex v during time period s 
and no detection at the vertex visited during time period s — 1. However, redehning 
the arc reward to depend not only on the arc’s head and tail nodes, but also on a 
previous node, ruins the longest-path structure of the bound-calculation problem. 
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A similar situation arises in vehicle routing problems for vehicles with 
turn-radius constraints or penalties. The classical approach to handle that situation 
is to duplicate each node a number of times equal to the number of nodes in the 
node’s reverse star. An arc in the resulting “node-expanded” network then carries 
information about three nodes, not only two, and a desirable network structure of 
the problem can be maintained. Fortunately, it is practical to carry out such a node¬ 
expansion approach in the problems of interest in this paper because the number of 
nodes in the reverse star is typically quite moderate. Hence, we proceed along the 
stated lines and develop a node-and-time expanded network, in which the improved 
static bound can be calculated by solving a longest-path problem. We refer to this 
improved bound as the directional static bound. 

For any n' G A/”, let lZ{n') d N he the reverse star of n', i.e., lZ{n') = 
{n G N'\{n,n') G A}. Then, for any n,n' G J\f\{h} such that (n,n') G A, we dehne 
an expanded node ^ = {n,n'). We do not expand the end node, so we set ^ = n. Let 
S be the set of all expanded nodes. Two expanded nodes G S are connected by 
an expanded arc if ■C = {n,n') and = {n',n"). Let the set of all expanded 

arcs be H. The node-and-time expanded graph is illustrated in Figure 4. 

We endow each expanded arc in the node-and-time expanded graph 
(S, H) with a reward similar to (11.10). To derive the exact form of this reward, we 
need the following building blocks. For any G V and t G T, let Mt{v,v') he a C 
by C identity matrix with the 0(n')-th diagonal element set equal to 1 — g{v,v',t)- 
We also let F(n') be the (;/)(n')-th column of the Markov transition matrix F. 

From (11.2) and the recursive application of (ILl), we see that the 
probability of detection along a path {vi}f^Q is given by 

q{{vi}f=o) = p{(j){vi), l)g{vo, ui, 1) + 
p(-, l)Mi{vo,vi)T{v2)g{vi,V2,2) + 
p{-, l)Mi{vo, vi)rM 2 {vi, V 2 )T{v 3 )g{v 2 , Vs, 3) -h (11.11) 
p(-, l)Mi(no, ni)FM2(ni, n2)FM3(n2, ^'3)^(n4)fl'(^'3, 4) -f- 
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Figure 4. A node-and-time expanded network from the time-expanded graph (Figure 
3). 

p{-, l)Mi{vo,Vi)TM2{Vi,V2) • . . . ■ TMT-l{vT-2,VT-l)'i^{vT)g{vT-l,VT,T). 

The expression (IF 11) gives insight into a class of bounds on the probability of detec¬ 
tion, including the static bound do{vt, t). If we replace Mi(-, ■) by the identity matrix 
in (II. 11), we hnd that 

q{{vi}f^o) < p(0(ni), l)^(no, Ui, 1) -f- 
p{-,l)T{v2)givi,V2,2) + 

p(-,l)FF(n3)^(n2,'y3,3) + (11.12) 

p{-, l)FFF(n4)5((n3,n4,4) -t- 

p(-,l)F^ ‘^T{vT)g{vT-i,VT,T). 
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In (11.12), the “reward” received during a time period is simply the expected number 
of detection during that time period and depends only on the current and previous 
vertices. Hence, it is possible to compute an upper bound on the optimal probability 
of detection by hnding a path {vi}f^Q that maximizes the right-hand side in (11.12). 
This calculation amounts to a longest-path problem and is, in fact, the approach 
in [34]. (Note, however, that [34] assumes that the glimpse detection probability is 
independent of the previous vertex.) 

If we replace each Mt{-, •) by the identity matrix everywhere except the 
last matrix of each line in (II. 11), we obtain 

q{{vi}f^o) < p{(j){vi), l)g{vo,Vi, 1) 

p{-, l)Mi{vo,vi)T{v2)g{vi,V2,2) + 

p(-, l)TM2{vi,V2)T{v3)g{v2,V3,3) + (11.13) 

p{-, l)TTM3{v2,V3)T{vi)g{v3,V4, A) + 


pi': l)r'^ ^Mr-i(nT-2, VT-i)T{vT)givT-i, vt, T). 

Now, the reward received during each time period also depends on the searcher’s 
position two time periods ago, and the problem of hnding a path that maximizes 
the right-hand side is no longer a longest-path problem. However, the bound remains 
valid with the following minor modihcation, where the maximization of a matrix with 
a single element different from zero or one is simply the maximization of that element; 

qi{vi}f^o) < p(0(wi), l)g{vo,Vi, 1) + 

p(-,l) ( niax Mi(n,ni) ) r(n 2 )^(t'i, ^ 2 , 2) + 

\v£-R{vi) J 

p(-,l)r( max M 2 (n,n 2 )) r(n3)5((n2,^^3,3)-h (11.14) 

pi': l)rr ( max M 3 (n,n 3 )) T{v4)g{v3,V4,4:) + 

\v£TI{v3) I 
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p(-,l)r'^^( max MT-iiv,VT-i)]T{vT)givT-i,VT,T). 

\v£-RivT-i) J 

After this modification, we see that the reward during each time period only depends 
on the current and previous vertices. Hence, again, it is possible to compute an upper 
bound on the optimal probability of detection by solving a longest-path problem. In 
fact, this is exactly the approach we described in Subsection B.2a and it can be shown 
that the reward in the longest-path problem see (II. 10), can be deduced from 
(11.14). Specihcally, when the current subpath in (II.10) is {uq}, we have for arc 
{n,n') = {{v, s), {v', s -f 1)) G ^ that 

Cn,n'= ( max M4n",n)) r(n')fl'(n,n',s-h 1). (11.15) 

\v"£TZ{v) j 

Using similar arguments, we dehne the directional static bound as fol¬ 
lows. Clearly, 

q{{vi}f=o) < l)g{vo, vi, 1) + 

p{-, l)Mi{vo, vi)T{v2)g{vi, V2, 2) -h 
p(-,l) ( max Mi{v,vi)]tM 2 {vi,V 2 )T{v 3 )g{v 2 , Vs,3) + (11.16) 

\ven{vi) J 

p(-,l)r( max M 2 (v,V 2 )]TMs(v 2 ,vs)T(vi)g(vs,Vi,^) + 
\v&l{v2) J 


p(-,l)r^ ^ max MT-2(v,VT-2)]'i^MT-l(vT-2,VT-l)T(vT)g(vT-l,VT,T). 

\v£-R{vt-2) j 

Hence, we can compute an upper bound on the optimal probability of detection by 
hnding a path maximizes the right-hand side of (11.16). This calculation 

amounts to a longest-path problem in the node-and-time expanded graph (2,12). 
The arc reward in this longest-path problem is deduced from (11.16). Specihcally, an 
expanded arc (^,^') = ({n",n), {n,n')) G H, with n" = {v",s — 1), n = {v,s), and 
n' = {v', s •+• 1), is endowed with the reward 


cw = p(-, i)r 


s-2 


max 


Ms_s(v"',v") VM,(v",v)V(v')g(v,v',s + 1). (11.17) 
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We refer to the node-and-time expanded graph (S, ^2) with the arc rewards from 
(11.17) as the node-and-time expanded network. Since the node-and-time expanded 
graph is acyclic, longest-path problems are solvable by standard shortest-path algo¬ 
rithms. 

In view of the above discussion, we obtain the following result. 

Proposition 1 For any n' G V and t E T, let 

(i) do{v',t) be the value of the longest-path from node {v',t) to node h in the 
time-expanded graph {Af^A) with are rewards given by (11.15), and 

(a) 6o{v,v',t) be the value of the longest-path from expanded node ((n,t —1), {v't)) 
to expanded node f in the node-and-time expanded graph (S, hi) with are re¬ 
wards given by (11.17). 

Then, both do{v',t) and So(v,v',t) are upper bounds on the probability of detection 
during time periods t -\- l,t -\- 2, ...,T for any path {vi}JLq with Vt-i = v and Vt = v'. 
Moreover, So(v,v',t) < do(v',t). 

We refer to So(v,v',t) as the directional static bound and see from 
Proposition 1 that it is at least as strong as the static bound. We demonstrate in an 
empirical study below that it may be substantially stronger. 

Clearly, building the node-and-time expanded graph (2,12), computing 
the associated rewards, and calculating the longest-paths take some computing time. 
However, the process is only carried out once before the start of Algorithm 1 and 
the computed bounds are stored for later use. Hence, the time for computing the 
directional static bounds remains small compared to the overall run time of Algorithm 
1. We conjecture that the use of the node-and-time expanded graph (2,12) with 
dynamic reoptimization of bounds will not be efficient due to the small but signihcant 
effort required to build the node-and-time expanded graph, compute associated arc 
rewards, and carry out the longest-path calculations. We observe that this conjecture 
appears to be aligned with [32], which alludes to the inefficiency of a dynamic bound 
based on more than one-time-step look-behind. 



In Table 2, we report computational results for Algorithm 1 with the 
directional static bound applied to the same problem instances as in Table 1. Columns 
3 and 4 give run times and number of branching attempts, respectively, for Algorithm 
1 using the directional static bound. We observe that, on average, the number of 
branching attempts has been reduced to 58.1% by using the directional static bound 
compared with the static bound. (Compare column 4 in Table 1 with column 4 in 
Table 2.) Similarly, the run times have been reduced to 58.2% by using the directional 
static bound. Since the reduction in branching attempts and run time are essentially 
identical, we conclude that the time for computing the directional static bound is 
small compared to the overall run time. 


g{v,v',t) 

P 

Algo. 1: 

Time 

(sec.) 

D-Static 

Branching 

attempts 

Algo. 2: 

Time 

(sec.) 

Static & Red. 
Branching 
attempts 

Algo. 2: 

Time 

(sec.) 

D-Static & Red. 
Branching 
attempts 


0.3 

2.59 

1,642,619 

0.22 

129,990 

0.19 

91,198 

0.3 

0.6 

1.80 

1,125,929 

0.17 

94,307 

0.16 

65,485 


0.9 

6.30 

4,032,951 

0.58 

354,677 

0.41 

223,435 


0.3 

3.50 

2,245,784 

0.31 

194,425 

0.27 

128,526 

0.6 

0.6 

1.45 

902,409 

0.17 

92,997 

0.14 

57,015 


0.9 

29.67 

19,321,387 

2.17 

1,442,612 

1.37 

844,747 


0.3 

3.59 

2,282,989 

0.34 

209,160 

0.28 

138,598 

0.9 

0.6 

1.66 

1,037,829 

0.17 

101,216 

0.17 

61,005 


0.9 

80.50 

52,527,302 

4.95 

3,323,101 

2.89 

1,837,647 


0.3 

3.77 

2,396,764 

0.36 

219,217 

0.28 

137,063 

0.99 

0.6 

1.72 

1,080,459 

0.19 

102,763 

0.17 

62,890 


0.9 

87.99 

57,427,410 

5.39 

3,592,696 

3.09 

1,974,871 


Table 2. Run time and number of branching attempts (counted in Step 4) for Al¬ 
gorithm 1 on problem instances of Table 1 using directional static bound (D-Static) 
and for Algorithm 2 using static bound and network reduction (Static & Red.) and 
directional static bound and network reduction (D-Static & Red.). 


d. Network Reduction 

In some practical situations, the searcher’s and the target’s initial po¬ 
sitions are relatively far from each other. Hence, for a number of time periods the 
searcher will only examine cells where the target is guaranteed not to be located. In 
these initial moves, the searcher is simply positioning itself for the later search. This 
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Time period 1 


Time period 2 




Figure 5. An area of interest composed of 5 by 5 cells. For each time period, the 
unshaded, lightly-shaded, heavily-shaded, and completely-shaded cells describe the 
regions where neither searcher nor target stay, only target possibly stays, only searcher 
possibly stays, and target and searcher possibly stay, respectively. 

is the situation in the numerical examples of [32]. In this subsection, we derive a 
network reduction technique that utilizes this situation. 

Consider the time-expanded graph Suppose that the searcher 

flies, for the hrst s time periods, along a subpath {nt}f=Q, with Ut = such 

that p{(l){vt),t) = 0 for all t < s, and p(0(t's),<s) > 0. Then, the searcher cannot 
detect the target prior to time period s. We refer to the last node Ug of the subpath 
{nt}f^Q as a node-of-£rst-contact. It is typically straightforward to determine a set of 
nodes-of-first-contact by applying a standard search algorithm ([1], pages 73-77) on 
the time-expanded network. In Figure 5, we illustrate a situation where the searcher 
can at the earliest detect the target during time period 4, and also note that the time 
period of all nodes-of-£rst-contact {vg, s) E J\f is s > 4. 
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We observe that there might be several different subpaths rit = 

{vtit) to reach a node-of-£rst-contact {vs,s) G J\f. However, the subpath used to 
reach {vs,s) is of no or little importance. In fact, if g{vs-i,Vs, s) does not vary 
with the choice of n^-i, all subpaths to {vs,s) have probability of detection equal 
to p{(l){vs), s)g{vs-i,Vs, s)] with p{(l>{vs),s) = p(-, l)r^“^r(ns). Thus, it is enough to 
find a subpath to each node-of-first-contact {vg, s) using a standard search algorithm, 
connect initial node {vq, 0) to {vg, s) directly with a “jump” arc representing the 
move to {vs,s), and ignore time periods l,2,...,s — 1 during the branch-and-bound 
algorithm. Clearly, this procedure may reduce the amount of branching attempts 
significantly. We can generalize this to the case with g{vg-i, Vg, s) varying with respect 
to ns_i, as discussed in Subsection C.2c. 

If a lower bound q on the optimal probability of detection is available in 
advance of the network reduction procedure described above, it may not be necessary 
to consider all nodes-of-first-contact. Fortunately in the single searcher problems 
(SSP), a lower bound is trivially obtained by computing the probability of detection 
along the path corresponding to the static bound (io('i^o,0). Furthermore, for each 
node-of-hrst-contact {vg,s), an upper bound on the probability of detection after 
time period s with no prior detections (e.g., the static bound do{vg,s)) is available 
in advance of both branch-and-bound and network-reduction procedures. The upper 
and lower bounds can be used to eliminate some nodes-of-first-contact that cannot 
lie on an optimal path. Specifically, if p{(l){vs), s)g{vg_i,Vg, s) -|- do{vg,s) < q, we can 
eliminate {vg, s) and all incoming and outgoing arcs. Thus, the amount of branching 
attempts may be reduced further. 

In order to take advantage of this reduced network in the branch-and- 
bound algorithm, we construct a second algorithm (Algorithm 2) that generalizes 
the branch-and-bound mechanism of Algorithm 1. Before we describe the algorithm, 
we need to clarify some notation. After the algorithm is presented, we describe the 
network reduction procedure in detail. 
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Given a subpath {nj}|hQ, ni = k E T, let /C(i), as before, be 

the set of triplets {vt,t,q{vt,t)) representing extensions of yet to be explored. 

Now, the index i of /C(i) refers to the depth from node {vo, 0) to the next node {vt,t) 
in the branch-and-bound tree. Since (wo^O) is directly connected to each node-of- 
hrst-contact {vs,s), the corresponding depth in the branch-and-bound tree is 1. For 
reporting purposes, a subpath is stored for each node-of-£rst-contact {vs,s). 

Algorithm 2 (Solves SSP). 

Input: A time-expanded network with arc length Cn^ or a node-and- 

time expanded network (S, G) with arc length The initial target 

distribution p{-, 1), Markov transition matrix F, and upper bound q. 

Output: An optimal search path V* = {vi}f^Q and its value q*. 

Step 0. Calculate So{v,v',t) for alH G T and v,v' eV such that {v,v') E S 
or do{v',t) for all t G T and v' E V, and calculate a lower bound q. Set 
i = 0, A(i) = {(uo,0, cx))}. 

Step 1. If )C{i) is empty, replace i by i — 1. Else, go to Step 3. 

Step 2. i = 0, stop: the last saved path is optimal and q is its probability 
of detection. Else, go to Step 1. 

Step 3. Remove from lC{i) the triplet (vt, f, f)) with the largest bound 

Step 4. If q{vt,t) < Q, go to Step 1. (Current subpath is fathomed.) 

Step 5. If i = 0, replace i by i + 1, and go to Step 3. (/C(l) is populated in 
the network reduction procedure.) 

Step 6. If t ^ T, then for each vertex v E .F('U^), calculate q(^vA T 1) from 
(II.5) using bounds do{v, f-|-l) or So(vt, + and add (u, f + 1, q(v,t + i)) 
to /C(i -|- 1). Replace i by i -|- 1, and go to Step 3. Else, let q = q{vt, t) and 
save the incumbent path {vi}JLq, and go to Step 1. 

We now present the network reduction procedure that can be imple¬ 
mented as part of Step 0 of Algorithm 2. The procedure assumes that a static bound 
do{v', t) is available as well as a lower bound on the optimal probability of detection q. 
If the directional static bound is available instead of the static bound, replace do{v', t) 
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by So{v,v',t) in the procedure below. We note again that the network reduction pro¬ 
cedure is valid under the assumption that g{vs-i, Vs, s) does not vary with the choice 
of Us-i £ T^{vs) among all nodes-of-£rst-contact {vs,s). We generalize this in Section 

C. 


Network Reduction Procedure Rl. 

Input: A time-expanded network (A/*, A) with arc length Cn,n', the initial tar¬ 
get distribution p{-, 1), Markov transition matrix F, and upper bound q. 

Output: Nodes-of-hrst-contacts {vs,s) which are not dominated by others. 

Step 1. Find all nodes-of-hrst-contact {vg, s) G A/”. If none exist with s > 1, 
then stop. 

Step 2. For each {vg, s), calculate q{vs, s) = p{(l>{vg), s)g{vg^i, Vg, s)+do{vg, s). 

Step 3. Eliminate all nodes-of-hrst-contact {vg, s) with q{vg, s) < q. 

Step 4. For each nodes-of-hrst-contact {vg, s) not eliminated, store the triplet 
{vg,s,q{vg,s)) in }C{1). 

Table 2 illustrates the ehect of the network-reduction technique as ap¬ 
plied to the same problem instances as in Table 1. Columns 5 and 6 of Table 2 
present the run time and number of branching attempts, respectively, for Algorithm 
2 with static bound and network reduction. On average, the run times and and the 
branching attempts are reduced to 5.3% and 5.0% of the corresponding numbers ob¬ 
tained without the network reduction technique (see columns 3 and 4 in Table 1), 
respectively. When applying both network reduction and directional static bound, 
we obtain the run times and numbers of branching attempts reported in columns 7 
and 8 of Table 2. It is clear that network reduction and directional static bound have 
complementary positive ehect and the run times and numbers of branching attempts 
are further reduced. 

Table 3 presents computational results for a larger problem instance 
with 15 by 15 cells and a time horizon T = 20. Again, the searcher starts in the 
upper-left cell and the targets starts in the center cell. As seen from Table 3, the run 
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times remain rather short for Algorithm 2 with directional static bound and network 
reduction (columns 3 and 4) while the times increases dramatically for Algorithm 1 
with dynamic bound (columns 5 and 6). Furthermore, Algorithm 2 with directional 
static bound and network reduction is less sensitive to the detrimental effect of a 
near stationary target (e.g., p = 0.9) and high glimpse detection probability (e.g., 
g{v,v',t) = 0.9 or 0.99), a case where the bounds tend to be weak. 


g{v,v',t) 

P 

Algo. 2: 

Time 

(sec.) 

D-Static & Red. 
Branching 
attempts 

Algo. 1: 
Time 
(sec.) 

Dynamic Bound 
Branching 
attempts 


0.3 

3.08 

103,811 

20.24 

328,672 

0.3 

0.6 

2.94 

66,185 

21.22 

311,645 


0.9 

3.58 

238,089 

107.11 

1,352,503 


0.3 

3.14 

122,941 

20.28 

248,727 

0.6 

0.6 

2.92 

51,977 

30.47 

288,738 


0.9 

4.41 

571,649 

701.17 

8,668,034 


0.3 

3.17 

129,276 

29.61 

292,818 

0.9 

0.6 

2.95 

57,353 

45.05 

404,299 


0.9 

5.94 

1,076,948 

2323.48 

30,173,994 


0.3 

3.13 

128,776 

31.97 

301,498 

0.99 

0.6 

2.92 

60,449 

48.42 

441,664 


0.9 

5.77 

1,016,918 

3092.63 

45,329,829 


Table 3. Run times and number of branching attempts for Algorithm 2, with direc¬ 
tional static bound and network reduction, compared with Algorithm 1 with dynamic 
bound on 15 by 15 cell search problem with time horizon T = 20. 


The same problem instances were also examined in [32], which reports 
a run time of 10.9 seconds for the case with g{v,v',t) = 0.6 and p = 0.6 using a 
C-|--|- implementation of Algorithm 1 with a dynamic bound running on a 2.6 GHz 
computer. We observe that our implementation appears to be somewhat slower than 
the one achieved in [32] with 30.47 seconds compared to 10.9 seconds (on a presumedly 
slightly slower computer). However, Algorithm 2 with directional static bound and 
network reduction appears to offer a noticeable advantage over Algorithm 1 with 
dynamic bound as derived in [32]. In principle. Algorithm 1 with dynamic bound can 
also be speeded up by using the proposed network reduction technique. However, we 
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have not examined that possibility. We adopt a directional static bound with network 
reduction as the basis for extension to the case with side constraints discussed in the 
next section. 

C. BRANCH-AND-BOUND ALGORITHM FOR RSSP 

We now turn the attention to the full problem with side constraints, i.e., the 
resource-constrained single searcher problem (RSSP) formulated in Section A. We 
first develop a static bound based on Lagrangian relaxation that can be used within 
a branch-and-bound algorithm in the form of Algorithms 1 and 2. Second, we briefly 
discuss the development of a Lagrangian directional static bound. Third, we develop 
a series of network reduction techniques. Fourth, we combine the resulting procedures 
and present the complete algorithm. 

1. Lagrangian Static Bound 

Consider the time-expanded network (A/”, A), see Subsection B.l, with the arc 
rewards Cny given in (11.15). Now, we also endow each arc (n, n') E A, n = {v,t — 1) 
and n' = with weights ri^n,n' = G X. While computing the static 

bound in the case of no side constraints amounts to solving a longest-path problem 
on the time-expanded graph, a similar bound in the case with side constraints will 
need to account for those constraints. Specihcally, a static bound can be obtained by 
solving a constrained longest-path problem. We formulate this problem as an integer 
program on a slightly modihed time-expanded graph. 

We use the same time-expanded graph as in Subsection B.l, with the following 
modihcations. We recall that V is the set of vertices where the search can end. Now, 
every node n = (v, t), w G V, t G T is connected to the artificial destination node h 
with an arc. This modihcation allows the searcher to terminate the search prior to 
time period T to avoid violating the side constraints. Furthermore, all arcs (n, n) with 
n = {v,T), V ^ V, are removed from the time-expanded network. This modification 
makes the searcher end at w G V. We still let A denote the set of all arcs. 
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We formulate the constrained longest-path problem on the time-expanded 
graph (A/*, A) as an integer program. We consider an ordering of A and let A denote 
the lA/”! by |Al| node-arc incidence matrix for the time-expanded graph. For each arc 
a = {n,n') G A, let = 1, A„',a = —1, and = 0 for any n" G Af,n" ^ n, n' 

be the elements of A. Let b denote the |A/’|-vector with = lAh = —1 and 6^ = 0 
for all n G A/'\{no,h}. We also dehne the additional notation; f = (ri,r 2 ,... i 
where T denotes the transpose. We collect the rewards Cn,n' in the |A|-dimensional 
row vector c. Also, for each i G X, we define the |A|-dimensional row vector r* to 
contain the weights ri^n,n' and we let R be the |X| by |A| matrix with as its rows. 
Finally, we let x be a |A|-dimensional column vector, where Xn,n' = 1 if arc {n,n') is 
used by a path, and zero otherwise. Then, the constrained longest-path problem, see 
[1], can be written as: 


= max cx 

(11,18) 


S.t. Ax = b 


Rx < r. 

(11.19) 


In principle, the solution of the constrained longest-path problem provides a static 
bound. However, the constrained longest-path problem is NP-complete even for the 
case with an acyclic graph ([22], pages 213-214). Hence, we prefer to avoid solving such 
problems within an algorithm. We proceed by introducing an additional relaxation 
that is motivated by the solution approach for the constrained shortest-path problem 
in [24, 9, 10]. 

Using the standard theory of Lagrangian relaxation (see, e.g., [1], Chapter 16) 
and an |X|-dimensional row vector A > 0, we hnd that 

< ^(A) = max cx — A(Rx — f) (11.20) 

xe{o,i}l-^ 

s.t. Ax = b. 
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Rewriting the objective function, we can optimize the Lagrangian upper bound 2 :(A) 
through 


2 ; = min 2 :(A) (11.21) 

A>0 

= min max (c —AR)x + Af (11.22) 

A>o xe{o,i}i-^ 

s.t. Ax = b. 

For any hxed A > 0, computing the upper bound ; 2 (A) simply requires the solution of 
a longest-path problem with Lagrangian-modihed arc lengths in an acyclic graph. The 
outer minimization over A can be solved by several methods [19, 5, 14, 9, 10]. Since 
we anticipate only a small number of side constraints, it suffices to use a repeated 
coordinate search. Given an optimal or near-optimal A, we obtain the Lagrangian 
static bound by carrying out one backward longest-path calculation in the time- 
expanded graph from node h to all nodes n & M using the Lagrangian modified 
arc reward c — AR. More specihcally, an arc (n,n') = ((v,s),(n',s -f 1)) G A, 
s = 1, 2,..., T — 1, is endowed with the reward 

Cn,n' ^ i,n,n' ; (11.23) 

i£l 

where Cn,n' is given by (11.15). 

We also derive and implement a Lagrangian directional static bound similar to 
the one in Subsection B.2c. However, the derivation is a straightforward combination 
of Subsection B.2c and the approach described above. Hence, we omit it. This 
derivation results in a similar Lagrangian problem to the one in (11.21), but now 
defined on the node-and-time expanded network. We still refer to the Lagrangian 
multiplier as A and the Lagrangian upper bound as 2 :(A). We denote the Lagrangian 
directional static bound computed in this way by 6o{v,v',t). Since the Lagrangian 
directional static bound is at least as strong as the Lagrangian static bound, we carry 
out the Lagrangian relaxation only in the node-and-time expanded network to hnd a 
Lagrangian multiplier A > 0. 
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2. Network Reduction 

We propose and examine three techniques for reducing the size of the network 
prior to application of a branch-and-bound procedure. First, we use dominance rules 
to eliminate edges that cannot be on an optimal path. Second, we describe the 
application of “preprocessing” techniques frequently used prior to solving constrained 
shortest-path problems. Third, we modify the procedure described in Subsection 
B.2d. 

a. Edge Dominance 

There are several situation where a vertex v' G Elv) can be eliminated 
as candidate for visit from vertex v. Such “dominance tests” are case dependent, but 
can be effective in reducing the number of edges. We describe one situation where we 
use “edge dominance.” 

In many practical situations, there are two resources: risk and fuel. 
If higher altitude implies lower risk and lower glimpse detection probability, and 
climbing to higher altitude consumes more fuel than level flight, then we can eliminate 
some edges in the graph (V,T) by edge dominance. Suppose that fi{v,v',t) and 
/ 2 (u, v', t) represent risk and fuel, respectively. Also suppose that the risk /i(u, u', t) = 
/i(u'), i.e., only depends on v'. Let 'ipiv) be the altitude of waypoint u G V. Then, if 
we have the above described situation, we use the following (one-step) procedure to 
reduce the size of the graph (V,T). 

Edge Dominance Procedure R2. 

Step 1. Delete any edge (u, v') G £ that satisfies /i(u') = 0 and 'ip{v) < ipiy'). 

We note that Procedure R2 takes advantage of the fact that if there is no risk at u', 
then there is no need to increase altitude when moving from v to v'. The altitude 
can be increased later if need be. 
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b. Preprocessing 

It is well known that the side constraints (11.19) can be used, prior 
to any main calculations, to identify nodes and arcs in the time-expanded network 
{N',A) that cannot lie on any feasible path [4, 15, 9, 10]. Such nodes and arcs 
can be eliminated from which reduces the size of the problem that needs to 

be considered when computing bounds and carrying out branching. Moreover, the 
reduction of the time-expanded network typically also strengthens the Lagrangian 
relaxation, i.e., reduces the gap z — z*, (see (11.21) and (11.20)), and, hence, reduces 
the need for branching. We adopt the follow procedure, adapted from [9, 10], to carry 
out arc preprocessing: 

Preprocessing Procedure R3. 

Input: A time-expanded network {J\f,A) and arc lengths 

Output: The updated (or reduced) time-expanded network {J\f,A). 

Step 1. Set number of iterations k and k = 1. 

Step 2. For alH G X and n G A/”, compute a minimum-weight rio-n subpath 
distance Riin) and a minimum-weight n-fi subpath distance rj(n) in (A/”, A) 
with respect to weights 

Step 3. Delete any arc (n, n’) G A with 

Ri{n) + + D(n') > fi for any i G T. (11.24) 

Step 4. If k < k and at least one arc was deleted in Step 3, replace khy k + 1, 
and go to Step 2. Else, stop. 

If a lower bound on the probability of detection is available, we also 
carry out similar preprocessing with respect to arc reward Cny and the Lagrangian 
modihed arc reward Cny = Cny — J2i£j K^i^ny, see [9, 10]. 

We describe the preprocessing procedure for the time-expanded net¬ 
work and argue that it improves the Lagrangian static bound. However, the same 
methodology applies to the node-and-time expanded network and it improves the 
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Lagrangian directional static bound. In our main algorithm (described in Subsection 
C.3), we also apply preprocessing to the node-and-time expanded network and denote 
that procedure R3'. 


c. Vertex Dominance for Distant Target 

As in Subsection B.2d, we consider the case where the searcher’s and 
the target’s locations are initially some distance apart and derive a network reduction 
procedure that takes advantage of that situation. In contrast to Subsection B.2d, 
the subpath used to reach a node-of-£rst-contact is now important since the resource 
consumption along different subpaths may be different. Hence, {vq, 0) now needs to be 
connected to each node-of-hrst-contact {vg, s) with multiple “jump” arcs representing 
the different possible subpaths and resource consumptions used to reach {vs,s). A 
standard path enumeration algorithm (see, e.g., [11]) can enumerate the different 
subpaths, at least as long as s is relatively small. Multiple arcs to a node-of-hrst- 
contact can also be used to model the situation with edge-dependent glimpse detection 
probability, a case ignored in Subsection B.2d. 

After all the arcs are generated to all the nodes-of-first-contact, a num¬ 
ber of them can be deemed uninteresting and can be eliminated using dominance rules 
of the form: If an arc from {vq, 0) to {vg, s) has no larger probability of detection and 
no smaller consumption of each resource as another parallel arc and the two arcs are 
not identical, the hrst arc is dominated and can be eliminated. In sets of identical 
parallel arcs, we also eliminate all but one. Trivially, arcs with resource consumption 
greater than the specified limits are also removed. Moreover, if a lower bound on 
the optimal probability of detection exists, it can be used to eliminate more arcs as 
described in the following. 

Below we formally describe this network reduction procedure based on 
vertex dominance for a distant target. We note that the procedure is more effective 
after (i) applying network reduction procedures R2, R3 and R3', (ii) hnding A that 
(approximately) optimizes the Lagrangian upper bound ; 2 (A) on the node-and-time 
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expanded network, and (iii) computing the directional static bound So{v,v',t) and 
the Lagrangian directional static bound So{v,v',t) for each node {v',t) G J\f. So 
we assume that these calculations have been carried out. During these calculations, 
feasible paths may be obtained. Such paths provide lower bounds on the optimal 
probability of detection. Let q denote the largest lower bound found so far. 

Vertex Dominance Procedure R4. 

Input: A time-expanded network {Af,A), arc lengths c„^„/ and Tn^n'- 

Output: The nodes-of-£rst-contact {vg, s) which are not dominated by others. 

Step 1. Find all nodes-of-hrst-contact {vs,s) G A/”. If none exist with s > 1, 
then stop. 

Step 2. For each node-of-first-contact {vs,s), enumerate all subpaths (no,0) 
to {vs,s). 

Step 3. For each node-of-first-contact {vs,s) and subpath V = carry 

out the tests: 

If any of the following is true, then eliminate V: 

(i) For some remaining {vo,0)-{vs, s) subpath V', ViiV) > rj('P') for all 
i G X and q{V) < q{V'). 

(ii) q{V) + 5o{vs-i, Vg, s) < q. 

(iii) r-iiV) + s)) > fi for some i G X. 

Step 3 can also be augmented with a test on the Lagrangian-modified probability of 
detection using 6o{v,v',t) if a near-optimal multiplier A is available. 

3. Algorithm 

We now state the complete algorithm for RSSP. The algorithm starts with 
network reductions procedures R2, R3, and R3'. The next step solves the Lagrangian 
problem (11.22) and determines a near-optimal A. (The calculations are actually 
carried out in the node-and-time expanded network as we prefer to use the Lagrangian 
directional static bound.) If a feasible path becomes available during the procedures 
described above, network reduction procedure R3' is repeated now using checks with 


41 



respect to arc reward and its Lagrangian modified arc reward. The next steps 
are to compute the directional static bound of Subsection B.2c (i.e., So{v,v' ,t)), the 
Lagrangian directional static bound as described in Subsection C.l (i.e., 6o{v,v' ,t)), 
and bounds on resource consumption along any path extension (i.e., rj(n), i G X). The 
final steps before the branch-and-bound procedure is to implement network reduction 
procedure R4. 

We implement the branch-and-bound procedure as an implicit path-enumeration 
in the time-expanded network. The procedure amounts to a depth-first search cou¬ 
pled with optimality and feasibility checks using the computed bounds, which follows 
[9, 10]. The complete algorithm takes the following form: 

Algorithm 3 (Solves RSSP). 

Input: A time-expanded network (A/*, A) with arc lengths Cn^n' and Vny and a 
node-and-time expanded network (S, hi) with arc lengths and The 
initial target distribution p(-,l), Markov transition matrix T, and upper 
bound q. 

Output: An optimal search path V* = {vi}f^Q and its value q*. 

Step 1. Apply network reduction procedures R2, R3 and R3'. 

Step 2. Find A that approximately optimizes the Lagrangian upper bound 
; 2 (A) in the node-and-time expanded graph (S, hi). If a feasible solution is 
found, set the probability of detection on the corresponding path equal to 
q. Otherwise, set q = —oo. 

Step 3. If a feasible solution is found so far, implement R3' also with respect 
to arc reward and Lagrangian modified arc reward using q. 

Step 4. Ignoring side constraints, compute the directional static bound 50 ( 1 ^, v', t) 
for all expanded nodes ^ = (n, n') in (S, O), with n = (n, t — 1), n' = {v\ t), 

V, v' e V. 

Step 5. Using A from Step 2, compute the Lagrangian directional static 
bound (5o(n,U,t) for all expanded nodes ^ = {n,n') in (S,r2), with n = 

{v, t — 1), n' = {v', t), V, v' e V. 

Step 6. For each i G X, compute the minimum distance rj(n) from each node 
n & M back to n by solving a single, backwards, shortest-path problem in 
the time-expanded graph {J\f,A) starting from h using arc lengths ri^n,n'- 
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Step 7. Apply network reduction procedure R4. 

Step 8. Apply a standard path-enumeration procedure (e.g., [11]) in {M^A) 
with the following modihcations: 

(i) The path-enumeration commences from no, but extends a current sub¬ 
path {niYi^Q along arc {nt,n) = ((nt,t),(n,t -|- 1)) if and only if the 
following conditions hold; 

• For all i G /, {no, ni,..., nt, n} can be extended to a path whose i-th 
resource does not exceed A, i.e., 

ri({no, ni,..., rit, n}) + rj(n) < A- (11.25) 

• (no, rii, can be extended to a path with probability of de¬ 

tection exceeding q, i.e.. 


g({no, ni,..., n*, n}) -h v,t + l) > q- (11.26) 

• (no, rii, can be extended to a path whose Lagrangian-modified 

probability is no less than q, i.e., 

t 

g({no,ni, ...,nt,n}) - 

iex 1=1 

- E + Ar 5o(^'i, V,t+Y>q- (11.27) 

i£l 

(ii) Whenever the algorithm identifies a path V with q{V) > q and ri{V) < 
fi,i G X, replace q by q(V). 


In Step 8, the checks (11.25), (11.26), and (11.27) prevent the enumeration 
of paths that can be determined, using the computed bounds, to not be optimal. 
Specihcally, (11.25) prevents the extension of subpaths that cannot result in a feasible 
path with respect to the side constraints. Since v,t + Y is a valid upper bound 
on the probability of detection during time periods f -|- 2, t -|- 3, ..., T, the left-hand 
side of (11.26) is an upper bound on the probability of detection along any path that 
starts with the subpath (no, ni,..., n^, n}. Hence, the subpath cannot be extended 
to a path with larger probability of detection than q if (11.26) fails. In (11.27), the 
probability of detection along {no,ni, is modihed by Lagrangian terms and 
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the Lagrangian directional static bound is used. The resulting check can be shown 
to be valid using standard Lagrangian relaxation theory and the argument above. 

We also implement Step 8 with a “branching strategy” based on the La¬ 
grangian directional static bound. Specihcally, we first consider extending the current 
subpath {riiYi^Q along the arc {nt,n) with the largest Lagrangian directional static 
bound among all the nodes in the forward star of rit. Second, we consider extend¬ 
ing with the node corresponding to the second largest Lagrangian directional 

static bound, etc. This branching strategy is analogous to the one in Step 3 of Algo¬ 
rithms 1 and 2. We have also experimented with using the directional static bound 
instead of the Lagrangian directional static bound and found it to usually result in 
comparable run times. However, the Lagrangian directional static bound appears 
faster, on average. 

Since Algorithm 3, in the worst case, enumerates all feasible paths, it is guar¬ 
anteed to hnd an optimal solution of RSSP. 

D. NUMERICAL EXAMPLE 

This section describes computational experiments with Algorithm 3 applied 
to RSSP with side constraints on risk exposure and fuel consumption. We carry out 
all experiments on the same computational platform as in Section B. 

We consider a military planning situation where a UAV is assigned a mission 
to search and detect a high-value moving target. Planners wish to determine a flight 
path over the area of interest (AOI) that maximizes the probability of detecting the 
target. The UAV will start its path at a known waypoint with a known fuel supply, 
and will return to the same waypoint before the fuel tank is empty. Doctrine specihes 
that the UAV cannot be assigned a path with higher risk than a specihc threshold. 
The AOI is partially under enemy control and any aircraft flying over the AOI could be 
shot down by enemy surface-to-air missiles (SAMs), anti-aircraft artillery, and small- 
arms hre. Flying at a high altitude would reduces that risk, but it will also reduce the 
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quality of UAV’s sensor readings. Consequently, the UAV may change altitude during 
the course of the mission to balance risk and sensor quality. Changing from low to 
high altitude consumes more fuel than level flight. Hence, the number of time periods 
available for search depends on the fuel consumption and therefore the vertical flight 
prohle. 

We model this situation by dividing the AOI into 10 by 10 cells (see Figure 
6 ). The airspace over each cell is vertically discretized into two altitudes (’Tow and 
high”). The heavily shaded cells (Ci) in Figure 6 represent an urban area over which 
the UAV’s risk is high and its glimpse detection probability is low. The unshaded cells 
(C 3 ) represent open terrain where there is no risk and the UAV’s glimpse detection 
probability is high. The lightly shaded cells (C 2 ) represent an area with intermediate 
risk and intermediate glimpse detection probability. We note that, while this problem 
instance does not directly relate to any real-world operational situation, we believe 
that it illustrate a fairly practical situation. 

We assume that the risks at different edges along a path are independent. If 
the probability of the UAV surviving edge (u, v') G ^ is a(v, v'), then the probability 
of surviving the path is simply nfLia(vi-i, vi). Let d be a lower limit on 

the survival probability. Then, a standard logarithmic transformation leads to the 
following constraint 

k 

- log o-(ui_i, Vi) < - log d, (11.28) 

i=i 

which is in the form (II.3) and (II.4) with fi{vi-i,vi,l) = —\oga{vi-i,vi) and A = 
- log d. 

For this computational experiment, we assume that the glimpse detection 
probability and the survival probability for an edge (u, v') G S depend only on the 
cell and altitude corresponding to vertex u' G V as listed in Table 4. We note that 
glimpse detection probability at high altitude is assumed to be 70% of that at low 
altitude and the failure probability (complement of the survival probability) at high 
altitude is 30% of that at low altitude. 
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Figure 6 . An area of interest composed of 10 by 10 cells and two altitudes. Heavily 
shaded cells (Ci), lightly shaded cells (C 2 ), and unshaded cells (C 3 ) describe risky, 
moderately risky, and non-risky area, respectively. The circle indicates the cell over 
which searcher starts, and the triangle specifies the initial position of the target. 

The UAV enters the airspace at high altitude over the northwest cell (cell 1; 
cells are numbered from left to right, and from top to bottom) and will return to 
the same cell at either high or low altitude at the end of the mission. The searcher 
is located at one vertex each time period and searches the corresponding cell. For 
the next time period, the searcher can stay at the same vertex, change altitude over 
the same cell, or move to a vertex (at any altitude) corresponding to a vertically or 
horizontally adjacent cell. The maximum number of time periods is T = 40, but the 
fuel-consumption constraint may limit the number of periods to less than 40. We 
assume that the fuel consumption at each time step is as follows: 10 units if there 
is no altitude change, 12(9) units if changing from low(high) altitude to high(low) 
altitude. The initial position of the target is the center of the high risk region (cell 
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Cells 

Altitude 

Glimpse probability 

Survival probability 

Cl 

low 

0.20 

0.960 


high 

0.14 

0.988 

C2 

low 

0.40 

0.980 


high 

0.28 

0.994 

Cs 

low 

0.60 

1.000 


high 

0.42 

1.000 


Table 4. Glimpse detection probability g{v,v',t) and survival probability (7{v,v') for 
different cells and altitude. 


68 ). The target remains in the current cell with a probability p = 0.6 for the next 
time period or moves to one of the vertically or horizontally adjacent cells with equal 
probability. 

The survival-probability limit is a threshold that is set by the commander or 
planner. A search path with lower survival probability than the threshold would not 
be accepted. In this experiment, we consider the survival probability limits 0.95, 
0.90,...,0.70, and fuel consumption limit 300, 325, ..., 450. 

We solve this problem instance using Algorithm 3. Tables 5, 6 and 7 report 
computational results for different combinations of survival probability and fuel limits 
for the UAV. When the fuel limit is tight (e.g., 300 and 325), the UAV cannot operate 
for the full duration of 40 time steps. We observe that increasing the fuel limit beyond 
425 does not increase the probability of detection as the time limit of 40 periods 
becomes active. The average run time is 580 seconds, with a standard deviation of 
792. All problem instances are solved within one hour and typically in much less. 
Figure 7 shows the optimal path given survival probability limit 0.90 and fuel limit 
400. The solid lines and the dashed lines represent flight segment at low and high 
altitude, respectively. 

We also consider the case with edge-dependent glimpse detection probability. 
Consider the same situation as earlier described, but now assume that a move to a new 
waypoint results in a lower glimpse detection probability than if the searcher already 
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was at that waypoint. Specifically, if n = v', we let the glimpse detection probability 
g{v,v',t) be as in Table 4; otherwise we replace g{v,v',t) by 0.1g{v,v',t). Figure 8 
shows an optimal path found given survival probability and fuel limits of 0.90 and 
400, respectively. In contrast to the case with edge-independent glimpse detection 
probability (Figure 7), the searcher now tends to stay for multiple time periods at the 
same waypoints in high-probability regions to reap the benehts of the corresponding 
high glimpse detection probability. The run times (not reported in detail) for the 
case with edge-dependent glimpse detection probabilities are, on average, 53 seconds, 
with a standard deviation of 71 seconds. The reduction in run time compared to 
the edge-independent case is caused by the often lower glimpse detection probability 
( 0 .l 5 f(n, F, t)), which tightens the bound. 


Fuel 

limit 

Survival prob. limit = 
Prob. Survival Fuel 

Detection Prob. 

0.95 

Run time 
(sec.) 

Fuel 

limit 

Survival prob. limit = 
Prob. Survival Fuel 

Detection Prob. 

0.90 

Run time 
(sec.) 

300 

0.131501 

0.962466 

300 

26.89 

300 

0.135098 

0.915157 

300 

29.58 

325 

0.153228 

0.952892 

321 

98.20 

325 

0.166933 

0.901925 

322 

32.66 

350 

0.176511 

0.952892 

350 

3313.64 

350 

0.194440 

0.915157 

350 

115.94 

375 

0.204686 

0.952892 

371 

661.43 

375 

0.217277 

0.901925 

372 

84.00 

400 

0.236651 

0.962466 

400 

664.96 

400 

0.246767 

0.902268 

400 

537.71 

425 

0.237081 

0.952892 

401 

2721.68 

425 

0.249996 

0.901925 

402 

361.60 

450 

0.237081 

0.952892 

401 

2724.22 

450 

0.249996 

0.901925 

402 

361.52 


Table 5. Computational results (probability of detection, survival probability, fuel 
consumption and run times) for Algorithm 3 with survival probability limit = 0.95 
and 0.90. 
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Fuel 

limit 

Survival prob. limit = 
Prob. Survival Fuel 

Detection Prob. 

0.85 

Run time 
(sec.) 

Fuel 

limit 

Survival prob. limit = 
Prob. Survival Fuel 

Detection Prob. 

0.80 

Run time 
(sec.) 

300 

0.145309 

0.851528 

300 

25.31 

300 

0.145809 

0.805953 

299 

22.03 

325 

0.173000 

0.851528 

320 

37.81 

325 

0.173218 

0.839535 

319 

37.39 

350 

0.203183 

0.851528 

350 

67.69 

350 

0.204891 

0.805953 

349 

64.17 

375 

0.222393 

0.851528 

370 

115.24 

375 

0.223377 

0.805953 

369 

116.69 

400 

0.255481 

0.851528 

400 

510.68 

400 

0.255813 

0.827710 

399 

880.65 

425 

0.255481 

0.851528 

400 

508.50 

425 

0.255813 

0.827710 

399 

880.40 

450 

0.255481 

0.851528 

400 

508.96 

450 

0.255813 

0.827710 

399 

880.46 


Table 6. Computational results (probability of detection, survival probability, fuel 
consumption and run times) for Algorithm 3 with survival probability limit = 0.85 
and 0.80. 


Fuel 

limit 

Survival prob. limit = 
Prob. Survival Fuel 

Detection Prob. 

0.75 

Run time 
(sec.) 

Fuel 

limit 

Survival prob. limit = 
Prob. Survival Fuel 

Detection Prob. 

0.70 

Run time 
(sec.) 

300 

0.150092 

0.753377 

300 

21.33 

300 

0.150277 

0.742767 

299 

21.38 

325 

0.175850 

0.753377 

320 

37.50 

325 

0.176068 

0.742767 

319 

31.16 

350 

0.205935 

0.753377 

350 

63.84 

350 

0.206200 

0.742767 

349 

53.42 

375 

0.223703 

0.753377 

370 

125.97 

375 

0.224003 

0.713056 

369 

121.33 

400 

0.255813 

0.827710 

399 

1220.30 

400 

0.255813 

0.827710 

399 

1280.52 

425 

0.255813 

0.827710 

399 

1217.21 

425 

0.255813 

0.827710 

399 

1277.52 

450 

0.255813 

0.827710 

399 

1218.91 

450 

0.255813 

0.827710 

399 

1282.42 


Table 7. Computational results 
consumption and run times) for 
and 0.70. 


(probability of detection, survival probability, fuel 
Algorithm 3 with survival probability limit = 0.75 
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Figure 7. An optimal path for survival probability limit 0.90 and fuel limit 400. The 
solid lines and the dashed lines represent flight segments at low and high altitude, 
respectively. (See Figure 6 for a description of the underlying figure.) 
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Figure 8. An optimal path for a case with edge-dependent glimpse probability, sur¬ 
vival probability limit 0.90, and fuel limit 400. The solid lines and the dashed lines 
represent flight segments at low and high altitude, respectively. (See Figure 6 for a 
description of the underlying figure.) 
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III. MULTIPLE-SEARCHER PROBLEM 

This chapter addresses the optimization of multiple searchers with path- and 
time-constraints (MSP). We aim to determine a path for each searcher that maximize 
the probability of detecting a moving target within a mission time by at least one 
searcher. We refer to such a set of paths as a “search plan.” We start by describing 
and formulating MSP. Then, we present three algorithms (an exact procedure and 
two heuristics) to hnd an optimal, near-optimal, or “good” search plan. The chapter 
also includes computational studies. 

A. PROBLEM DESCRIPTION AND FORMULATION 

The multiple-searcher problem (MSP) is an extension of SP which considers 
a hnite set of searchers J = {1,2, ...,J}. In this section, for ease of reference, we 
formulate MSP by retracing the relevant portions of Section A in Chapter II. 

The area of interest is discretized into a hnite set of cells C = (1,..., C*} and the 
time horizon is discretized into a hnite set of time periods T = {1,2, ...,T}. A target 
occupies one cell each time period and moves among cells according to a Markov 
process with known transition matrix P. Let p(-,t) = [p(l, t),p(2, t),... ,p(C, t)], 
where p(c, t) is the probability that the target is in cell c G C at the beginning of time 
period t G T and the target has not been detected before t by any searcher. We refer 
to p(-,t) as the undetected target distribution. The initial target distribution p(-, 1) 
is known. 

The multiple searchers move through a designated airspace over the area of 
interest with the goal of hnding the moving target on the ground. The airspace over 
each cell c G C is vertically discretized into a set of altitudes Td = {1,2,..., H}. For 
any c G C and h eH, we refer to the cell-altitude pair (c, h) as a waypoint where the 
searchers can loiter and carry out search of cell c. We model the designated airspace 
by a directed network (V,T), with set of vertices V and set of directed edges S, in 
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which vertices v = {c, h) E V represent waypoints and directed edges e = {v, v') G £ 
represent transition between waypoints n, v' G V. The searchers can only transit 
between two waypoints that are located physical adjacent to each other. Let J^{v) dV 
be the set of vertices that are adjacent to n G V. We adopt the convention that 
V G J^{v) for all n G V. Then, the set of edges £ = {(n,n') G V x V| n' G 

During each time period t E T, each searcher is located at a particular vertex 
(waypoint). Any number of searchers can occupy a vertex in the same time period. 
We also assume that there is no transit time between waypoints. Hence, (n, v') E £ 
simply represents search at waypoint v followed by search at waypoint v' in the next 
time period. We note that the edge (n, v) E £ represents searching at waypoint v for 
two consecutive time periods. 

Let 0 ; V —> C be the function that specifies the cell over which a vertex 
is located, i.e., cell 0(v) is searched from vertex v. We denote the initial vertex 
(waypoint) of searcher j E J prior to time period Ihy Vq E V. We also define W C V 
to be a set of possible destination vertices for searcher j E J . When we describe an 
arbitrary searcher, we may omit the superscript and simply write vq and V. 

For any k E T and Vi E V, I = 0,1,2,k, such that {vi-i,vi) E £ for 

all I = l,2,...,k, let the sequence V = denote a directed Vo-Vk subpath. If 

Vk E V, then the directed VQ-Vk subpath is a directed VQ-Vk path. In this notation, 
a searcher flies from vq to some Vk E V along a directed vo-Vk path. For a specific 

searcher j E J, we denote its directed v^-Vk (sub)path as Thus a 

search plan for the fleet of searchers is described asV = {V^,V^, 

We adopt the following target-detection model. If the target is in cell c G 
C during time period t E T and only searcher j E J is at the same time at a 
waypoint above cell c, i.e., = c, then detection occurs with a glimpse detection 

probability g^{vl_i,vl,t), where vl_i G V is the waypoint of searcher j during time 
period t — 1. Hence, the glimpse detection probability during time period t depends 
on the previous and current waypoints for the searcher. We assume that the glimpse 
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detection probability of each searcher is mutually independent. In this notation, the 
probability of detecting the target in cell c G C by any searcher during time period 
t E T and no prior detections becomes 

Q{c,t) = p{c,t) 1- Yl {1- , (III.l) 

j&J\(l>{vl)=c 

We refer to Q{-,t) = [Q{l,t),Q{2,t),... ,Q{C,t)] as the reward vector at time period 
t. 

Since p{-,t) is the undetected target distribution at the beginning of time 
period t, it depends on searches up to time period t — 1. Specifically, if p{-, t) is given 
and each searcher searches a cell from its waypoint during time period t E T, the 
undetected target distribution at the beginning of the next time period t + 1 is 

p(.,t + l) = {p{;t)-Q{;t)]T. (III.2) 

Thus, the probability of detection for search plan V is defined as 

( 111 . 3 ) 

t=l c=l 

We refer to the problem of maximizing (III.3) as the multiple-searcher problem (MSP). 

B. ALGORITHMS FOR MSP 

We develop one exact algorithm and two heuristics to solve MSP. The exact 
algorithm is a straightforward extension of the branch-and-bound (B/B) algorithm 
in Chapter II. However, the B/B algorithm for solving MSP requires an expanded 
network structure to account for multiple searchers. Thus, in each time period t E T, 
we consider the combined location Vt = (n/, n/,..., n/) of all searchers. We call 
this combined location a configuration. We treat a configuration in the same way 
we treated a vertex in SP, and construct a time-expanded (configuration) network. 
The B/B algorithm for MSP is implemented in this time-expanded network. In that 
network, a search plan is simply a sequence of conhgurations {Vi}f^Q, k E T. We also 
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use the notation where ni = Thus we refer to a search plan also as 

a “conhguration path.” When the meaning is clear from the context, we refer to a 
conhguration path as a path. 

We present two heuristics to solve MSP. The hrst one is a variant of the 
expected detection heuristic (HI) in [13]. In HI, a conhguration path (search plan) 
is generated by iteratively extending the current conhguration subpath by the hrst 
arc in the longest-path, with respect to the expected number of detections, from the 
current conhguration to the last time period. We note that a longest-path calculation 
is required every time the current conhguration subpath is extended. Hence, HI is 
similar to Algorithm 1 with dynamic bounds (see Chapter H) as they both require 
repeated longest-path calculations. As an alternative approach, we present a similar 
heuristic corresponding to the static bound (see Chapter H). We denote our approach 
the static bound heuristic (SBH). 

The second heuristic algorithm is using the Cross-Entropy (CE) method [38, 
39]. The CE method involves an iterative procedure where each iteration composes 
of two phases. First, the method generates random samples (i.e., search plans in 
our case) according to a probability distribution. Second, the method updates the 
parameters in the probability distribution based on a subset of the “best” samples, 
the so-called “elite” samples. This process increases the chances of an optimal or near- 
optimal solution to appear within the next set of samples. This dissertation appears 
to be the hrst one studying the CE method in the context of search problems. 

We can assume without loss of generality that each searcher’s path consists 
of T -|- 1 vertices. For simplicity of notations, we also assume that: (1) there is no 
end-point restriction for any searcher j G JT”, i.e., W = V; (2) the glimpse detection 
probability g^{v,v',t), j E J' is independent of v (i.e., g^{v,v',t) = g^{v',t)); (3) the 
airspace has only one altitude, i.e., there is only one vertex corresponding to each cell 
in the area of interest. It is straightforward to generalize the proposed algorithms to 
account for situations without these assumptions. 
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1. Branch-and-Bound Algorithm 

In Chapter II, we presented a specialized branch-and-bound (B/B) procedure 
for solving SB. We utilize the same enhancement schemes to develop a B/B algorithm 
for MSP. However, the B/B algorithm for solving MSP requires an expanded network 
structure to account for multiple searchers. We consider a directed graph (Vc, 
where Vc is the set of possible conhgurations and Ec is the set of directed edges 
representing possible transitions in one time period between two conhgurations. Thus 
a search plan can be described as V = {VjiLo where (V/-i,Vz) ^ Ec for all I = 
1,2, ...,T. We note that a search plan is also described (in the previous section) as 
V = {V\ ...Pb where = K}f=o- 

The description of the B/B algorithm for solving MSP follows the presentation 
of the algorithms for SP presented in Chapter II. Given a conhguration subpath 
{ViYiZq, t G T, let /C(t) be the set of triplets of the form q{Vt,t)) representing 
extensions of {V;};=o to be explored. The hrst element V) refers to the next 
conhguration to form, the second element t is the time period to form the conhguration 
Vt, and the third element q(yt,t) is an upper bound on the probability of detection 
along any conhguration path that starts with the conhguration subpath {V)}f=o- The 
upper bound g(V),t) consists of three parts. Let be an upper bound on 

the probability of detection during time periods t + l,t + 2, ...,T, given that the 
conhguration of time t is Vt, and no detection occurs along the conhguration subpath 
The two other parts are the probability of detection on the conhguration 
subpath {Vi}\Zq and the probability of detection during t. Hence, 

q{Vu t) = qiiVYzY + J2Qivt) + d,(V„ t). (HI.4) 

cec 

We also let q denote the largest detection probability found so far among all the ex¬ 
amined conhguration paths. 
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Consider a configuration subpath {V;}*^q, t eT, and let Pg{-,t) be the unde¬ 
tected target distribution after search along i.e., 

Pgi-,^) = (III.5) 

For any integer s > t, s E T, we also define 

Pr{;s;t)=pg{;t)r-K (IIL6) 

As seen, pr{c,s;t) is the probability that the target is in cell c at time period s > t 
and there was no detection during search along the subpath Hence, target 

distribution pr{-,s]t) at time period s > t ignores the effect of search after time 
period t. If the subpath is {Vq}, i.e., t = 0, we dehne for notational convenience 

pr(-,s;0)=p(-,l)r-\ (III.7) 

for any s > 0, s G T. Moreover, we dehne pr{c, t;t) = 0 for all c G C and t = 0,1,..., T. 

We construct a time-expanded conhguration graph (A/'c, Ac) from the graph 
(Vc, Sc) in the same manner as the development of the time-expanded network for 
SP in Chapter II (for details, we refer to Section B in Chapter II). For any integer 
k < T +1 and nodes ni = {Vi,l) E Me, I = 0, l,...,/c, such that (nj_i,nj) G Ac for 
all I = we let the sequence {ni}f^Q denote a conhguration subpath in the 

time-expanded conhguration graph {Me, Ac). 

For some t E {0,1,...,T — 1}, suppose that a conhguration subpath 
in the graph (Vc, Sc) is given. We endow each arc (n, n') = ((V, s), {V, s -1-1)) G Ac, 
s = t,t + 1,...,T — 1, in the time-expanded conhguration graph {Me, Ac) with a 
“reward” 



- 

- 




Pr{d,s + l;t) - 

Qids)T{c,d) 


n - 9^s + ^) 

c'GC 

- 

ce7^(c') 




(III.8) 


where F(c, d) is the c-d element of the Markov transition matrix F. Thus this time- 
expanded conhguration network {Me, Ac) has the same structure as the time-expanded 


58 



network (A/*, A) of Chapter II and the longest paths in this network provides the static 
bound dolVt, t). The B/B for MSP is implemented in this time-expanded conhguration 
network {Nc,Ac). 

We consider generalizations of directional static bound and network reduction 
from Chapter II. To obtain a directional static bounds, we would need to generate a 
node-and-time expanded “conhguration” network corresponding to the node-and-time 
expanded network for SP. Since the number of arcs in the node-and-time expanded 
conhguration network becomes extremely large (e.g., a problem instance with 5 by 5 
cells, 3 searchers, and time horizon 7 has more than one hundred million arcs), we 
realize that the directional static bound technique is not computationally practical for 
MSP. Hence, we adopt the static bound in our B/B algorithm for MSP. The difficulty 
of network size occurs also in the time-expanded conhguration network when we 
consider large-size problems, which means our B/B algorithm is practical only for 
small-size problems. 

It is straightforward to generalize the network-reduction technique based on 
the initial positions of the searchers and the target of Chapter II to MSP. The required 
“conhgurations”-of-hrst-contact (14, s) G A// are now dehned as the hrst conhgura¬ 
tion where at least one searcher can obtain contact with the target. Thus our B/B 
algorithm utilizes this network reduction technique. The resulting algorithm takes 
the following form: 

Algorithm 4 (Solves MSP). 

Input: A time-expanded conhguration network (A//, Ac) with nodes n G A/” 
and arcs (n,n') G A where n = {V,t — 1), n' = Arc lengths c„^„/ 

in (A//, Ac), the initial target distribution p{-, 1), Markov transition matrix 
P, and upper bound q. 

Output: An optimal conhguration path V* = {V)}^o value q*. 

Step 0. Calculate the static bound do(y,t) for all t G T and V G Vc, and 
calculate a lower bound q. Set i = 0,/C(i) = {(Vq; 0, cx))}, where Vq is 
the initial conhguration. Implement network reduction procedure R5 (see 
below). 
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Step 1. If /C(f) is empty, replace ihy i — 1. Else, go to Step 3. 

Step 2. \i i = 0, stop: the last saved configuration-path is optimal and q is 
its probability of detection. Else, go to Step 1. 

Step 3. Remove from lC{i) the triplet (yt,t,q(yt,t)) with the largest bound 

Step 4. If qiVtit) < q, go to Step 1. (Current subpath is fathomed.) 

Step 5. If f = 0, replace f by f -1- 1, and go to Step 3. (/C(l) is populated in 
the network reduction procedure.) 

Step 6. If t R, then for each conhgnration R G .^(V)), calculate q(V^t 1) 
from (III.4) using bounds do{V,t + 1) and add {V,t + l,q{V,t + 1)) to 
/C(f -|- 1). Replace f by f -|- 1, and go to Step 3. Else, let q = qiVt^t) and 
save the incumbent conhgnration path {V/}^g, and go to Step 1. 

Network Reduction Procedure R5. 

Input: A time-expanded conhgnration network {Me, Ac), the initial target 
distribution p(-, 1), Markov transition matrix E, and upper bound q. 

Output: The conhgurations-of-hrst-contacts (R,s) which are not dominated 
by others. 

Step 1. Find all conhgurations-of-hrst-contact (V^,s) G Me- If none exist 
with s > 1, then stop. 

Step 2. For each (14, s), calculate g(I4, s) = I]cec Q(c, -s) -|- do{Vs, s). 

Step 3. Eliminate all conhgurations-of-hrst-contact (14, s) with g(I4, s) < q. 

Step 4. For each conhgurations-of-hrst-contact (14, s) not eliminated, store 
the triplet ( 14 , s, qiVs, s)) in /C(l). 

We implement Algorithm 4 and examine its performance on small-scale in¬ 
stances of MSP. The problem instances are as follows: An area of interest (AOI) 
consists of 5 by 5 cells. The number of searchers is 1, 2, and 3. The searchers, which 
have a constant glimpse detection probability 0.6, operate in the AOI and their moves 
are restricted to vertically or horizontally adjacent cells. The target remains in the 
current cell with probability 0.6 or moves to one of the vertically or horizontally ad¬ 
jacent cells with probability 0.4. The diherent moves are equally likely. All searchers 
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depart the upper-left corner cell, and the target starts at the center of the AOI. The 
time horizon is from 5 to 10. We carry out all experiments on the same computation 
platform as in Chapter II. 

Table 8 reports the run times for different combinations of number of searchers 
and time horizons. Column 2, 4 and 6 of Table 8 show the run times of Algorithm 
4 for the case of 1, 2 and 3 searchers, respectively. For the case with 3 searchers, 
the run time increases exponentially with increasing time horizon, and the case with 
time horizon 10 cannot be solved in several days. The case with 3 searchers and 
the time horizon 9 takes about 5 days to guarantee an optimal solution, however, 
the optimal solution is found after only 174.78 seconds. The remaining time is spent 
proving optimality by fathoming other possible conhguration paths. Table 9 presents 
the optimal probability of detection obtained by Algorithm 4. Column 2, 5 and 8 of 
Table 9 show the optimal probability of detection for the case of 1, 2 and 3 searchers, 
respectively. Figure 9 illustrate an optimal search plan for the case with 3 searchers 
and time horizon 9. 


Time 

Horizon 

1 searcher 
B/B SBH 
(sec.) (sec.) 

2 searchers 
B/B SBH 
(sec.) (sec.) 

3 searchers 

B/B SBH 

(sec.) (sec.) 

5 

0.00 

0.00 

0.02 

0.00 

4.17 

1.84 

6 

0.00 

0.00 

0.03 

0.02 

5.67 

2.17 

7 

0.00 

0.00 

0.06 

0.02 

58.89 

2.72 

8 

0.00 

0.00 

0.44 

0.02 

4,341.51 

2.91 

9 

0.00 

0.00 

5.19 

0.03 

426,874.92 

3.25 

10 

0.00 

0.00 

75.77 

0.03 

N/A 

3.70 


Table 8. Run times for the branch-and-bound algorithm (B/B) and the static bound 
heuristic (SBH) on 5 by 5 cell search problem with 1-3 searchers and time horizon 


T = 5-10. 


In view of Tables 8 and 9, it is clear that Algorithm 4 is practical only for 
small-size problems (with few searchers, short time horizon, and/or small area of 
interest). As noted earlier, the time-expanded configuration network (necessary for 
implementing the B/B procedure) cannot be generated for large-size problems since 
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Time 

Horizon 

1 searcher 
B/B SBH 

Rel. 

gap 

2 searchers 
B/B SBH 

Rel. 

gap 

3 searchers 
B/B SBH 

Rel. 

gap 

5 

0.306483 

0.306483 

0.00 

0.474213 

0.474213 

0.00 

0.579710 

0.579710 

0.00 

6 

0.351647 

0.351241 

0.00 

0.535954 

0.521669 

0.03 

0.643001 

0.622074 

0.03 

7 

0.389043 

0.380220 

0.02 

0.581175 

0.561550 

0.03 

0.691865 

0.679234 

0.02 

8 

0.416987 

0.404325 

0.03 

0.618416 

0.574542 

0.07 

0.728375 

0.711876 

0.02 

9 

0.444506 

0.426829 

0.04 

0.647400 

0.620582 

0.04 

0.754400 

0.739376 

0.02 

10 

0.465594 

0.438671 

0.06 

0.673168 

0.648007 

0.04 

N/A 

0.762183 

N/A 


Table 9. Probability of detection obtained by the branch-and-bound algorithm (B/B) 
and the static bound heuristic (SBH) on 5 by 5 cell search problem with 1-3 searchers 
and time horizons T = 5-10, and relative gap between SBH and B/B solutions. 


the extremely large number of data (arcs) can not be stored. Thus, to solve large-size 
problems, alternative heuristic approaches may be necessary. 

2. Static Bound Heuristic 

Dell et ah [13] present seven algorithms (one exact procedure and six heuris¬ 
tics) to solve MSP. Among their heuristics, the expected detection heuristic (HI) 
performs best under a broad range of conditions. As described in the beginning of 
this section, HI constructs a search plan (conhguration path) by extending a config¬ 
uration subpath one arc at the time. Each extension is determined by a longest-path 
calculation on a time-expanded conhguration network where arc lengths are given by 
the expected number of detections. 

We present an alternative approach, the static bound heuristic (SBH), which is 
similar to HI in [13] but requires only one longest-path calculation and uses improved 
arc lengths as described below. In Algorithm 4 (i.e., the B/B procedure for MSP), we 
calculate the static bounds in advance of the main calculations. The bound calculation 
corresponds to a longest path in the time-expanded conhguration network (A//, Me) 
with arc length given by (HI.8). That longest path, which specihes a search plan, 
is the solution provided by SBH. We observe that the longest path calculations in 
HI amounts to hnding a conhguration path with the largest expected number of 
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Figure 9. An optimal search plan for the 5 by 5 cell search problem with 3 searchers 
and the time horizon 9. 

detections. In contrast, the longest path calculation in SBH uses arc lengths (III.8) 
which effectively amounts to hnding the conhguration path with the largest expected 
number of detections while at each node along the path accounting for the effect of 
search at the previous node. As demonstrated in [32] and discussed in Section 1 of 
Chapter II, accounting for the previous node signihcantly improves the conhguration 
path found by the longest path calculation. 

Column 3, 5 and 7 of Table 8 report the run times for the cases of 1, 2 and 
3 searchers, respectively. Since SBH corresponds to a longest-path problem in an 
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acyclic network, the run times tend to be extremely short. Column 3, 6 and 9 of 
Table 9 show the obtained probability of detection, and column 4, 7 and 10 describe 
the relative gap to the corresponding optimal probability of detection. For each case, 
the obtained detection probability is quite close to the corresponding optimal value. 
Thus SBH performs well for these problem instances. However the performance tends 
to be worse as the problem size becomes large. Moreover, this SBH is not available 
for large-size problem since the time-expanded conhguration cannot be generated. 

3. Cross-Entropy Method 

The cross-entropy (CE) method was developed by Rubinstein [38, 39] for solv¬ 
ing rare event simulations and combinatorial optimization problems. The CE method 
derives its name from the cross-entropy (or Kullback-Leibler) distance between two 
probability distributions (e.g, an optimal importance sampling distribution and an 
estimated distribution). The CE method is an iterative procedure consisting of two 
steps: (1) random samples (i.e., search plans for our case) are generated according to 
a parameterized probability distribution, and (2) the generated search plans are eval¬ 
uated using the objective function (i.e., detection probability), and the parameters 
of the sampling distribution are updated based on the ’’elite” samples in a manner 
which increases the possibility that an optimal or near-optimal solution appears in 
the next iteration. 

The CE method is a global search heuristic, and is somewhat similar to ge¬ 
netic algorithms. However the CE method has a simpler scheme for the change of 
population generation parameter compared to genetic algorithms. Boer et ah [36] 
compare in details the CE method and other heuristics such as simulated annealing, 
genetic algorithm, and ant colony method. The CE method has been applied to a 
large number of combinatorial optimization problems [39, 40, 33, 44, 36, 12] but this 
dissertation appears to be the first one applying the CE method to search problems. 
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a. CE Algorithms 

In the multiple-search problem (MSP), a possible search plan is de¬ 
scribes as a sequence of conhgurations V = {Vz}^q where {Vi-i^Vi) G £c for all 
I = Thus, a crude way to hud a search plan is as follows. Let EiV) = 

{V G Vc| (fo, V) G £c} be the forward star of conhguration V. Start at the given ini¬ 

tial conhguration Vq = {vq,Vq, ... ,Vq). Select an arbitrary conhguration from JP(Vo) 
and denote it fo. Next, choose an arbitrary conhguration from EiVi) and denote it 
V 2 . The same process is repeated until the conhguration Vt is obtained. 

The CE method takes similar steps, but selects the next conhguration 
according to a probability distribution. For each node n G Me in the time-expanded 
conhguration network (A/),,Me), we dehne a probability distribution an,n' over the 
outgoing arcs (n,n') G Ac (i.e., Y,{n,n')^Ac^n,n' = !)• Then, given a conhguration 

subpath k E T, the conhguration subpath is extended by randomly selecting 

an arc {nk,n') G Ac according to the probability distribution <Jnk,n'- 

Clearly, if an,n' is degenerate at each node n G M^ i.e., there exists an n' 
such that (n, n') G Ac and an,n' = 1, then the probability distribution uniquely speci- 
hes a conhguration path. The CE methods aims to iteratively update the probability 
distribution so that it converges to a degenerate distribution specifying the optimal 
conhguration path of MSP. We formalize this approach in the next algorithm. 


Algorithm 5 (Basic CE Algorithm). 

Parameters. Sample size M, elite size stopping parameter s (e.g., 

s = 5), and smoothing parameter /3. 


Step 0. Calculate static bounds do(n) for all n E Me and a lower bound q. 
Set = q and i = 1. 


Step 1. For all n E Me and n' E A’(n), dehne probability distribution 
such that 

(i) _ doin') 


^£2in'(iT(n) doin ) 


(III.9) 


Step 2. Generate M search plans based on probability distribution a 


(i) 

n,n' * 
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Step 3. Choose elite samples among the generated search plans, and 

calculate the detection probability of the best elite If > g, then 
g = g*^*) 

Step 4. If g*^*) = gh“i) = ... = gh-^)^ then stop. Else go to Step 5. 

Step 5. Update probability distribution using the elite samples 

and parameter f3 as described below. Replace i by i + 1 and go to Step 2. 

In Step 0, the static bounds d^in) are calculated for all node n G A/'c in 
the time-expanded conhguration network. At this time, we obtain a lower bound g 
of the optimal detection probability (i.e., the detection probability along the conhg¬ 
uration path corresponding to do{no)). The best detection probability initially (0th 
iteration), denoted is simply g. In Step 1, we dehne the initial probability dis¬ 
tribution using the static bounds. Next, in Step 2, based on the probability 
distribution cr^n,., we generate M search plans Pi, 7^2, • • •, Pm randomly. After that, 
in Step 3, the M plans are evaluated using the objective function (i.e., detection prob¬ 
ability), and the best performing elite samples are extracted. If the detection 

probability of the current best elite sample g*^^^ is better than the current lower bound 
g (the best detection probability so far), g is updated. Step 4 stops the algorithm if 
the best values found in the previous several iterations were identical. We hope that 
the repetition of values is due to a near degenerate probability distributions and 
that the best found search plan is close to the optimal one. However, that is not guar¬ 
anteed. If the stopping criterion is not satished in Step 4, probability distributions 
are updated using the current elite samples as follows. Let be the number 

of times one of the elite samples use the arc (n, n') in the current iteration. The 
probability distribution is updated using both the current probability distributions 
and contribution of the elite samples by setting 

At'’ = + (1 - (iii.io) 

where 0.4 < {3 < 0.9 as suggested in [36] based on empirical studies. 
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We note that the parameters (M, s, j3) must be specihed to 

implement Algorithm 5. Boer et ah [36] suggest that the sample size M is chosen 
according to the size of the problem (e.g., M = rm, where m is the number of arcs in 
the time-expanded conhguration network, and r is a constant), and the elite sample 
size is taken to be approximately O.OIM. Boer et ah [36] also report that an adaptive 
choice of the parameter settings speeds up the convergence. Thus, we implement nu¬ 
merical tests for a variety of cases and develop an adaptive CE algorithm. Specifically, 
we set the elite sample size and smoothing parameter as = O.OIM and (3 = 0.4, 

respectively. Furthermore, we vary the sample size M in the range between 
and adaptively. M™'*"' is set according to the size of the problem instance, and 

j^max _ 

The CE method has asymptotic convergence properties. Specifically, 
under the assumption that the optimal solution is unique, Costa et ah [3] give neces¬ 
sary and sufficient conditions under which the optimal solution is approached, with 
probability one, as the number of iterations goes to inhnity. For MSP, the CE method 
provides good solutions (as shown below in numerical tests). However it does not 
guarantee an optimal solution since MSP often has multiple optimal solutions (e.g., if 
two searchers are identical) and since the CE method is implemented with a heuristic 
stopping criterion. The multi-extremal property of the MSP solutions may provide 
the following undesirable condition. In each iteration, probability distributions are 
updated using the elite samples including the best elite. Thus, in the next iteration, 
at least, the best elite of the previous iteration is expected to appear in the samples 
if the sample size is large. However, the multi-extremal property may bring unstable 
probability distributions and the current best elite may be worse than that in the 
previous iteration. In this case. Algorithm 5 is unlikely to stop in reasonable time. 
For that undesirable situation, we intentionally stop the algorithm and provide the 
best solution, but label it as “unreliable.” 
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In view of the above discussion, we develop a specialized adaptive CE 
algorithm to solve the MSP. In this algorithm, we refer to the sample size at 
iteration as Mj. 


Algorithm 6 (Adaptive CE Algorithm). 

Parameters. Minimum and maximum sample sizes and and 

smoothing parameter 13. 


Step 0. Calculate static bounds do{n) for all n G Afc and lower bound q. Set 
= q and i = 1. 


Step 1. For all n G Me and n' G JP(n), dehne probability distribution a, 


(i) 

n,n' 


such that 


= 

^n,n' 


doin') 


(IIMl) 


^2in'(iT(n) d^ijl ) 

Set Mi = 

Step 2. Generate Mi search plans based on probability distribution 


Step 3. Choose = O.OlMj elite samples among the generated search 

plans, and calculate the detection probability of the best elite If > 
g, then g = g^*^ 

Step 4. If gh) > gh-P^ update using the elite samples. Replace 

ihy i + 1. Set Mi = M™*” and go to Step 2. Else go to Step 5. 


Step 5. If the best detection probabilities are identical for three iterations 
(i.e., g^*^ = gh“^) = then stop. (The obtained solution is considered 

reliable). Else go to Step 6. 

Step 6. If Mj < increase the sample size Mi = Mi + 

and go to Step 2. Else go to Step 7. 

Step 7. If Mk = Mi_i = ... = Mi_^, then stop. (The obtained solution 
is considered unreliable). Else update using the latest elite 

samples. Replace ihy i + 1. Mi = M™'*"' and go to Step 2 


Algorithm 6 is based on the time-expanded configuration network, 
which is problematic for large-scale problems. However, we overcome this difficulty 
by considering a time-expanded network for each searcher. For each searcher j G 77, 
we construct a time-expanded network as in Chapter II. For each we dehne the 


68 



probability distributions We then generate a search plan by generating a path 

= {n/}^Q for each searcher j according to the probability distributions The 
search plan is simply V = We update the probability distributions 

j G J', by considering the joint search effect of all searchers. Specihcally, we 
determine as before the elite search plan and dehne be the number of times 

one of the elite samples use the arc (n, n') for search j in the current iteration. 
The probability distribution for each searcher j is then updated by setting 


_ a 

^n,n' — P 


^n,n' 


/i(b 

2^n'£T(n) ^n,n' 


+ (1 ~ /^)t 


J'h) 

n,n'' 


(III.12) 


An advantage of this approach is that we can treat large-scale problems since we 
avoid generating the time-expanded conhguration network. We summarize this “de¬ 
composition” approach next. 


Algorithm 7 (Adaptive CE Algorithm with Decomposition). 

Parameters. Minimum and maximum sample sizes and and 

smoothing parameter /?. 

Step 0. Generate the time-expanded network , j G J. In each , calculate 
static bounds d^iji) for all nodes n G . Generate a search plan V by 
collecting each searcher’s path corresponding to the static bound (io(no). 
Galculate detection probability q{V) from (III.3). Set q = q^^'> = q{V) and 
Z = 1. 

Step 1. For j G JT”, define the probability distributions by 


Set Mi = 


^ doin') 
En'er(n)diin'y 


(III.13) 


Step 2. Gonstruct Mi search plans by generating each searcher’s path based 
on probability distribution For each search plan, V calculate the 

detection probability qiV) from (III.3). 

Step 3. Ghoose = O.OlMj elite search plans among the generated 

search plans, and set the detection probability of the best elite as qd\ 
If gb) > then q = gb) 
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Step 4. If > update Vj G J using searcher j’s path con¬ 
tributing to elite samples. Replace f by i -f 1. Mj = M™'*"' and go 

to Step 2. Else go to Step 5. 

Step 5. If the best detection probabilities are same for three iterations (i.e., 
gW = g(*-i) = then stop. (The obtained solution is considered 

reliable). Else go to Step 6. 

Step 6. If Mj < increase the sample size 

and go to Step 2. Else go to Step 7. 

Step 7. If Mfc = Mj_i = ... = Mj_ 5 , then stop. (The obtained solution is 
considered unreliable). Else update j G 77, using (III. 12). Replace 

f by i -I- 1. Set Mi = M™*” and go to Step 2 

In Algorithm 7, for each j G 77, the probability distributions <J^n% 
are rather poor in the initial iterations compared to the situation in Algorithm 6. 
However, the distributions are gradually improved as the cooperative search effort is 
accounted for when determining the elite search plans. 

The decomposition approach is practical for large-size problems since it 
only requires the time-expanded network for each single searcher. Thus, the approach 
is easily be implemented using the node-and-time expanded network with the poten¬ 
tial for more effective generation of paths (see Chapter II for a discussion about the 
advantage of the note-and-time expanded network over the time-expanded network). 
The size of the node-and-time expanded network is larger than time-expanded net¬ 
work, but it is substantially smaller than the time-expanded conhguration network. 
Since the algorithm using the node-and-time expanded network is essentially identical 
to Algorithm 7 (except for using the node-and-time expanded network), we omit to 
explicitly describe that algorithm but refer to it as Algorithm 7N. 

h. Numerical Tests 

We implement Algorithm 6, Algorithm 7, and Algorithm 7N using the 
twelve problem instances listed in Table 10. The problem instances are similar to the 
ones in Tables 8 and 9, and they have the same assnmptions and parameter settings 
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(e.g., start positions, glimpse detection probability and target stationary probability, 
etc). The minimum sample size M™*" is decided according to the size of problem 
instance, see columns 5 and 6 of Table 10 where m is the number of arcs in the 
network corresponding to the algorithms. As before, we set = 5M™*" and 

(3 = 0.4. Since Algorithm 6 uses the time-expanded conhguration network, it is 
available only for small problem instances, i.e., cases A2, A3, B2, and B3. The CE 
heuristics (Algorithms 6, 7, and 7N) are compared with Algorithm 4 (i.e., the B/B 
algorithm). When available, we compare the optimal value from Algorithm 4 with the 
search plan obtained by the CE heuristics. We also compare with incumbent solution 
available in Algorithm 4 at the point in time when the CE heuristics terminate. 


Case 

Area 

size 

Time 

horizon 

Number of 
searchers 

Min. sample size 
Algo. 6 Algo. 7/7N 

A2 

5 by 5 

9 

2 

O.lm 

10000m 

A3 

5 by 5 

9 

3 

0.01m 

10000m 

B2 

5 by 5 

18 

2 

O.lm 

10000m 

B3 

5 by 5 

18 

3 

0.01m 

10000m 

C2 

15 by 15 

18 

2 

N/A 

1000m 

C3 

15 by 15 

18 

3 

N/A 

1000m 

D1 

15 by 15 

27 

1 

N/A 

100m 

D2 

15 by 15 

27 

2 

N/A 

100m 

D3 

15 by 15 

27 

3 

N/A 

100m 

El 

15 by 15 

28 

1 

N/A 

100m 

FI 

15 by 15 

29 

1 

N/A 

100m 

G1 

15 by 15 

30 

1 

N/A 

100m 


Table 10. Test cases with Algorithm 6, Algorithm 7 and Algorithm 7N with area size, 
time horizon, number of searchers, and minimum sample size, m is the number of 
arcs in the network used in the corresponding algorithm. 
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(a) Single searcher 

Tables 11 and 12 report the numerical results for problem instances with 
a single searcher. Columns 2 and 3 in Table 11 present the detection probabilities 
obtained by the CE heuristics and column 4 shows the corresponding optimal solution 
obtained by Algorithm 4. We hnd that Algorithms 7 and 7N have comparable solution 
quality. The last column in Table 11 reports the probability of detection of the 
incumbent search plan of Algorithm 4 at the time Algorithm 7 terminates. For 
example, for case Dl, the best detection probability available to Algorithm 4 at time 
50.98 seconds (see Table 12) is 0.304046. Table 12 reports the corresponding run 
times (columns 2 and 3 for Algorithms 7 and 7N, respectively, and column 5 for 
Algorithm 4). In Table 12, we have also included the time at which Algorithm 4 Ends 
the optimal solution (but not yet proven optimal), see column 4. The run time of 
Algorithm 4 is clearly slower than that of the CE heuristics. We also observe that the 
run time of Algorithm 7 is much faster than that of Algorithm 7N. In a comparison 
of columns 2, 3, and 5 in Table 11, we find that the CE heuristics terminates with 
a better search plan than what is available from Algorithm 4 at the same time. The 
solution quality of the CE heuristics is quite good for all the instances examined. 


Case 

CE heuristic 

Algo. 7 Algo. 7N 

Branch-and-bound (Algo. 4) 
Optimal Early term. 

Dl 

0.305254 

0.305254 

0.305254 

0.304046 

El 

0.311225 

0.313101 

0.313101 

0.307411 

FI 

0.320277 

0.320716 

0.320719 

0.313043 

G1 

0.325968 

0.325968 

0.327823 

0.320131 


Table 11. Probability of detection for problem instances with a single searcher. 


(h) Two searchers 

Tables 13 and 14 describe the computational results for problem in¬ 
stances with two searchers. The columns show the same content as the corresponding 
columns in Tables 11 and 12. Algorithm 4 performs well for the smallest problem 
instance (A2), but is not available for large instances. Algorithms 7 and 7N obtain 
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CE heuristic 

Branch-and-bound (Algo. 4) 


Algo. 7 

Algo. 7N 

Found optimal 

Proved optimal 

Case 

(sec.) 

(sec.) 

(sec.) 

(sec.) 

D1 

50.98 

226.08 

151.14 

752.32 

El 

47.28 

229.17 

536.40 

2732.36 

FI 

79.26 

302.88 

1955.64 

9690.79 

G1 

60.17 

309.08 

6255.47 

34820.26 


Table 12. Run times on problem instances with a single searcher. 


solutions in practical times and their solution quality is better than that of Algorithm 
6. In addition, they can obtain solutions for large problems. Hence, it appears that 
Algorithms 7 and 7N are superior to Algorithm 6 for problem instances with two 
searchers. Algorithm 7 and Algorithm 7N are comparable, but the former takes less 
times to obtain a solution. The run time of the CE heuristics is smaller for case D2 
than for C2, see Table 14. This is counterintuitive as cases D2 have a longer time 
horizon than C2. However, the randomness in the algorithms may cause such effects. 
Overall, it appears that Algorithm 7 is the algorithm of choice as in the case of two 
searchers. 


Case 

Algo. 6 

CE heuristic 
Algo. 7 

Algo. 7N 

Branch-and-bound (Algo. 4) 
Optimal Early term. 

A2 

0.646586 

0.647400 

0.647400 

0.647400 

0.647400 

B2 

0.799565 

0.801566 

0.801552 

N/A 

N/A 

C2 

N/A 

0.336465 

0.336483 

N/A 

N/A 

D2 

N/A 

0.474782 

0.476186 

N/A 

N/A 


Table 13. Probability of detection for problem instances with two searchers. 

(c) Three searehers 

Tables 15 and 16 report the numerical results for the cases with three 
searchers. The columns are the same as the corresponding columns in Tables 11 and 
12. We observe that Algorithm 7 dominates Algorithm 6 in solution quality and run 
times. Furthermore, it is available for the large-scale problem instances. As in the 
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Case 

Algo. 6 
(sec.) 

CE heuristic 

Algo. 7 Algo. 7N 
(sec.) (sec.) 

Branch-and-bound (Algo. 4) 
Found optimal Proved optimal 
(sec.) 

A2 

6.78 

11.09 

48.45 

0.67 

5.19 

B2 

35.23 

169.89 

517.05 

N/A 

N/A 

C2 

N/A 

650.93 

1666.43 

N/A 

N/A 

D2 

N/A 

210.43 

951.36 

N/A 

N/A 


Table 14. Run times on problem instances with two searchers. 


case of two searchers, Algorithms 7 and 7N are comparable, but the former tends to 
require less computing time to obtain a reasonable solution. 

For the smallest size problem instance (A3), Algorithm 7 obtains an 
optimal solution (however it is not guaranteed) in 49 seconds. On the other hand. 
Algorithm 4 obtains an optimal solution in 174.78 seconds and proves it optimal 
in about 5 days. The best detection probability after 49 seconds in Algorithm 4 is 
0.749631. Thus the CE heuristic Algorithm 7 appears to generate good solutions 
quicker than the B/B based Algorithm 4. 

Comparing the CE heuristic (Tables 11-16) with the Static Bound 
Heuristic (SBH) (Tables 8 and 9 second to last row), we hnd that the SBH dom¬ 
inates the CE heuristic in run time (3.25 seconds compared to 49 seconds for case 
A3). However, the SBH is available only for small-size problem. In addition, the so¬ 
lution quality of the CE heuristic is better than that of the SBH based on the limited 
tests . 

In view of the above results. Algorithm 7 appears to be an overall 
efficient heuristic for solving MSP. 
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Case 

Algo. 6 

CE heuristic 
Algo. 7 

Algo. 7N 

Branch-and-bound (Algo. 4) 
Optimal Early term. 

A3 

0.753688 

0.754400 

0.754400 

0.754400 

0.749631 

B3 

0.889145 

0.891720 

0.893104 

N/A 

N/A 

C3 

N/A 

0.436528 

0.436528 

N/A 

N/A 

D3 

N/A 

0.587159 

0.593178 

N/A 

N/A 


Table 15. Probability of detection for problem instances with three searchers. 


Case 

CE heuristic 

Algo. 6 Algo. 7 Algo. 7N 
(sec.) (sec.) (sec.) 

Branch-and-bound (Algo. 4) 
Found optimal Proved optimal 
(sec.) 

A3 

144.75 

49.00 

249.42 

174.78 

426,874.92 

B3 

767.39 

271.43 

1567.62 

N/A 

N/A 

C3 

N/A 

674.85 

3815.69 

N/A 

N/A 

D3 

N/A 

297.30 

1454.73 

N/A 

N/A 


Table 16. Run times on problem instances with three searchers. 
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IV. MULTIPLE HOMOGENEOUS 
SEARCHER PROBLEM 

This chapter focuses on the multiple-searcher problems where all searchers are 
identical, i.e., the multiple homogenous searcher problem (MHSP). The goal for the 
searchers is to minimize the probability of not detecting the target within a mission 
time, which is equivalent to maximize the probability of detecting the target within 
the same time. We utilize the convexity of the nonlinear objective function (the non¬ 
detection probability), and introduce an exact algorithm using cutting planes (outer 
approximations). Under certain assumptions, the problem becomes equivalent to a 
large-scale linear mixed-integer program. We also present several new cuts for MHSP 
and demonstrate their effect in computational tests. 

A. PROBLEM DESCRIPTION AND FORMULATION 

Chapter III presented three algorithms (the branch-and-bound procedure and 
two heuristics) to solve the multiple searcher problem (MSP), which are also applica¬ 
ble for solving MHSP. Here we introduce an alternative approach especially tailored 
to solve MHSP. This section formulates MHSP by following [7, 18]. We use the same 
notations and assumptions as in MSP. Specihcally, we assume that: (1) there is no 
end-point restriction for any searcher; (2) the airspace has only one altitude, i.e., 
there is only one vertex corresponding to each cell in the discretized area of interest; 
and (3) the search effect at the current waypoint is independent of the previous way- 
point. In addition, we assume that the “search effect” is described by an exponential 
detection function instead of an arbitrary function as in earlier chapters. That is, 
when the number of searchers in cell c at time period t is yc,t, the probability of de¬ 
tecting the target in that cell during that time period given the target occupies cell c 
at time period t is described as 1 — exp(—Q;c,ti/c,t) instead of 1 — (1 — g{c, f))^'"’* as used 
in Chapter HI. Here, the term ac^t is the detection rate for a single searcher in cell c 
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and time period t defined by = —log(l — g{c, t)). We assume that g{c, t) < 1 (i.e., 
the sensor is not perfect) so that ac,t is finite. Additionally, we let ui{t) be the target’s 
cell at time period t E T, the sequence of cells u = (cc;(l), ci;(2),..., cc;(T)) denote a 
target path, and be the probability that the target takes that path u. The set of 
all possible target paths is denoted as hi. In this notation, MHSP is formulated as 
follows. 

Formulation of MHSP 


Indices 

c, c' cells (c, c' G C = {1,..., C}) 
t time step (f G T = {0,1,..., T}) 

u target path {u G hi) 

Parameters 

ac,t detection rate for a single searcher in cell c and time period t 

ac,t if cell c is on target path u at time period t, zero otherwise. 
i/cfl number of searchers in cell c at time period 0 
Puj probability that the target takes path u 

Variables 

Xc^c't number of searchers that is redistributed from cell c in time 
period t to cell c' in time period t + I 
yc,t number of searchers in cell c in time period t 

Formulation 


min f\y) = 

E 

p^ exp - 

' E ^c,tyc,t 



V 

c,t 

s.t. 





-1 = 

y ) ^c,c' 

,t Vf = l,, 

c'£-R{c) 


c'eJ*^(c) 


yc,t = Y. 

^c,c^ 

,t VfGT 



c'eJ*='(c) 

Xc,c',t, Vcp integer 


We refer to this problem which minimizes the probability that the target is 
not detected during time horizon T subject to network flow-balance constraints as 
the multiple homogeneous searcher problem (MHSP). The objective function of this 
nonlinear mixed-integer program is clearly convex. For reference later, elements of 
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the gradient Vf{y) of the objective function f{y) are 

P<^ac,texp(- 2^ac,tl/c,t). 
Ol/c,t c,t 


(IV.l) 


Since the formulation uses all possible target paths cu G hi, if the number 
of possible target paths is extremely large, calculating the objective function value 
/(y) or its gradient Vf{y) for any solution y becomes extremely computationally 
expensive. Brown [7] introduces an alternative formulation for the case of Markovian 
target movement using conditioning on the target’s position at a time period as 
follows. 


Given a solution y, let Vc^tiy) be the probability of the target visiting cell c at 
time period t without being detected in time periods 1, 2, ..., t — 1, and be the 

probability that the target departs cell c in time period t and is not detected by any 
searches in time periods f + 1, t + 2, ..., T. Let r.^tiy) = Vi,t{y)iT 2 ,t{y)i ■ ■ ■ )Tc,t{y)] 
and s.^t{y) = S 2 ,t{y), • • •, sc,t{y)]- We dehne r.^i{y) = p{-, 1), where p(-, 1) is a 

given initial target distribution, and Sc^xiy) = 1 for any cell c G C. Let L be a target 
transition matrix. Thus r.^t{y) and s.^t{y) are calculated recursively by 


r.,t{y) = [ri,t-i{y) exp{-ai^t-iyi,t-i), • • •, rc,t-i{y) exp{-ac,t-iyc,t-iW, (IV.2) 

■S;tiy) = [si,t+iiy) exp(-ai,t+i|/i,t+i),..., rc,t+iiy) exp{-ac,t+iyc,t+iW\ (IV.3) 

where L' is the transpose matrix of L. Since we only consider a target moving ac¬ 
cording to a Markov transition matrix, we can apply this result. In this notation, for 
any t = 1, 2,..., T, the objective function is alternatively expressed as 


fiy) = IlG,t(l/)exp(-ac,tl/c,t)sc,i(l/), 

cec 


and elements of the gradient V f{y) are 

dfiy) 


dyc,t 


= -ac,trc,tiy) exp(-ac,t|/c,t)sc,t(|/)- 


(IV.4) 


(IV.5) 


Thus, for any solution y, the objective function value and the gradient can be evalu¬ 
ated with a moderate computational effort. 
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B. ALGORITHMS FOR MHSP 

As mentioned, MHSP is a convex nonlinear mixed-integer program. This sec¬ 
tion introduces an exact outer approximation (OA) algorithm to solve MHSP using 
the convexity of the objective function [30]. The OA algorithm refers to the fact that 
the surface described by a convex function lies above the supporting hyperplane to 
the convex function at any point. A supporting hyperplane at a point takes the 
form /(yW) + V/(|/(*))'(|/ - yW) (see Figure 10). The algorithm iteratively generates 
supporting hyperplanes (cuts) and accumulates them to provide successively improv¬ 
ing linear approximations of the nonlinear convex function (see Figure 10). The linear 
approximation underestimates the objective function. We start by describing a basic 
OA algorithm. After that we introduce several new cuts and present computational 
results in the subsequent subsections. 

1. Basic OA Algorithm 

The basic OA algorithm is described as Algorithm 8 (see below), which solves 
the following master problem as part of its calculations. Specihcally, we denote that 
the master problem of the /c-th iteration of Algorithm 8 for MPl(/c) and its optimal 
value and optimal solution for and respectively. 

Formulation of Master problem : MPl(/c) 


min 

s.t. 

Z > - y^^) i = 0,1, k-l 

c^G7^(c) 

yc,t ^ ^ ^ ^ 

c' eJ^(c) 

Xc,c',t, yc,t- integer 

Algorithm 8 (Basic OA Algorithm). 

Step 0. Set a relative optimality tolerance 5 > 0. Choose an initial feasible 
solution y^^\ 
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Figure 10. Linear supporting functions. 

Step 1. Calculate and (see (IV.4) and (IV.5)). Set lower 

bound q = 0, upper bound q = f{y^^^), and k = 1. 

Step 2. If the gap {q—q)/q < S, stop. Else, solve the master problem MP1(A;), 
and obtain its optimal value and optimal solution y^^\ If zi^) > q, then 
q = z^^\ 

Step 3. Calculate and Vf{y^^^) (see (IV.4) and (IV.5)). 

Step 4. If f{y^’^^) < g, then q = f{y^^^). Replace /c by /c + 1, and go to Step 

2 . 

As seen, Algorithm 8 generates one cut in each iteration (a strategy we refer 
to as “single-cut”) at a solution determined by the master problem MP1(A;). 

It is well known that using a nonzero optimality tolerance when solving, at 
least the initial master problems in Algorithm 8 may reduce overall computing time 
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as cuts can be generated more cheaply. We examine this possibility. Let > 0 be 
a relative optimal tolerance for the solution of the master problem MP1(A;). Then, 

(k) 

the value obtained from the master problem may not be a valid lower bound on 
the optimal value of MHSP. However, (1 — is a valid lower bound and we use 

that instead of in step 2 of Algorithm 8. We note that when is constant at 
zero for iterations after some hnite iteration. Algorithm 8 is an exact algorithm and 
provides an optimal solution after a hnite number of iterations [16]. 

We examine Algorithm 8 using the problem instance “case C3” in Table 10 
of Chapter 111 which involves 15 by 15 cells, a time horizon of 18, and 3 searchers. 
We also adopt the same assumptions and parameter settings as in Chapter 111. Our 
program is coded using the General Algebraic Modeling System (GAMS) Distribution 
22.5 [21] and is implemented on the same computational platform as in Chapters 11 
and 111. The master problem MPl(/c) is solved by the CPLEX 10.0 solver [27]. In the 
master problem, we assume that the optimality tolerance is constant with value 
0 (version 1) or with value 0.03 (version 2). 

All searchers depart the upper-left corner cell. Thus, as a choice of initial 
feasible solution (search plan), we consider the situation that all searchers loiter at 
the depot during the whole time horizon. We refer to this search plan as a “trivial 
search” plan. This trivial plan is clearly unwise, but we apply this plan as an initial 
feasible solution in the preliminary tests. Later we introduce a procedure to obtain a 
better initial feasible solution. 

In the initial numerical tests, the goal is to compare various approaches thus 
we terminate the calculations after one hour. Table 17 reports lower/upper bounds 
on the optimal value and the relative gap between the bounds at every ten minute 
during one hour calculations using Algorithm 8. As reference, the best known non¬ 
detection probability is 0.563472 (=1-0.436528) which is found in Table 15 (case C3) 
in Chapter 111. The one hour run time includes the solution time for the master 
problems (Step 2) and the overhead time for generating cuts, etc. (Steps 1 and 3). 
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Algorithm 8 versions 1 and 2 use 96.4% and 85.1% of the run time to solve the master 
problems, respectively. During one hour, version 1 executes only 25 iterations because 
the tighter optimality tolerance in the master problems, while version 2 executes 96 
iterations. Even though the cuts in version 2 may not be as “deep” as the one 
generated in version 1, it is clear from Table 17 that the extra effort to solve the 
master problem optimally may not be worth it. 

In these test, the relative optimal tolerance of the master problem is con¬ 
stant (0 or 0.03) for the one hour calculations. We note that if is not zero. 
Algorithm 8 may regeneration a cut. This can be prevented in many ways, especially 
if we access the master problem solver. However, in these tests, as well as in tests 
below, we adopt the simple approach of adjusting the relative optimality tolerance 
g(*:) Por example, e^^^=0.03 appears to be a practical number if the goal is a 5% 
near-optimal solution. 


Time 

Algo. 8 

version 1: e 

(fc) = 0 

Algo. 8 version 2: 

= 0.03 

(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.410705 

0.620493 

0.511 

0.467784 

0.576832 

0.233 

20 

0.422249 

0.599865 

0.421 

0.471399 

0.576832 

0.224 

30 

0.428358 

0.599865 

0.400 

0.473591 

0.576832 

0.218 

40 

0.437659 

0.585815 

0.339 

0.475205 

0.574215 

0.208 

50 

0.448581 

0.585815 

0.306 

0.475316 

0.574215 

0.208 

60 

0.449400 

0.585815 

0.304 

0.475316 

0.572823 

0.205 


Table 17. Numerical results for Algorithm 8 versions 1 and 2 for every ten minutes 
of calculations. The best known non-detection probability is 0.563472 (=1-0.436528) 
found in Table 15 (case C3) in Chapter III. 


2. New Cuts for OA Algorithm 

This subsection introduces several new cuts and demonstrates their effect. 
Initially two new cuts (multi-cut and strong-cut) are presented. After that we prove 
that, under certain assumption, MHSP is equivalent to a large-scale linear mixed- 
integer program, which motivates an approach for obtaining a good initial solution. 


83 



Next, two additional cuts (relative-cut and symmetric-cut) are presented. Finally, we 
develop a specific OA algorithm by combining these cuts effectively. 

(a) Multi-cut 

In Table 17, we learned that at least a moderate number of cuts is necessary 
to bring up the lower bound. The aim of the multi-cut presented here is that, in 
each iteration, we generate multiple cuts instead of single cut to accumulate many 
cuts fast and push up the lower bound quickly. The multi-cut utilizes a similar 
technique to the multicut version of the L-shaped method in stochastic programming 
with two-stage recourse (see Chapter 5 in [ 6 ]). The basic idea is that, for the objective 
function f{y) = PujfM where f^{y) = exp(-Ec,t we consider outer 

approximation of fujiy) for all a; G hi instead of for /(y), and generate |r 2 | cuts at 
each iteration. However, the number of possible target paths |r2| is extremely large. 
Thus, according to the target movement during the initial few time steps, we define 
U (typically U « |r2|) “scenarios” Hi, H 2 , •••, H^/, where a scenario is a subset 
of target paths (i.e., C H, m = 1,2,..., 7/). Moreover, each scenario is mutually 
exclusive (i.e., fl = 0 , m 7 ^ u') and each possible target path oj belongs to 
some scenario (i.e., = H). For example, from a known target location at time 

period t = 1 there are five possible target movements (stay in the initial position or 
go up/down/left/right) for the next time period. Thus we can define five {U = 5) 
scenarios based on the target movement conditioned during the initial two time steps. 
Similarly, for the initial three time steps (t=l,2 and 3), twenty five scenario {U = 25) 
are defined if boundary effects are ignored. Let be the probability that the target 
takes any path oj such that u & VLu. 

In order to apply this multi-cut approach, the objective function (IV.4) and 
the gradient (IV.5) need to be expressed differently. When the solution is y, let r^t{y) 
be the probability of target visiting cell c at time period t without being detected 
in time periods 1 , 2 , ..., t — 1 given the target takes a path oj corresponding to 
scenario u. Similarly, s^^{y) is defined as the probability that the target departs cell 
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c in time period t and is not detected by any searches in time periods t + 1, t + 2, 
T given the target chooses a path uj corresponding to scenario u. Let r^t{y) = 
..., and s'f^^{y) = [sltiy), ^Itiy), • • •, (?/)]• % considering 

the conditioning on the target’s position with respect to scenario u, both r^t{y) and 
s^^{y) are recursively calculated similar as r.^t{y) and s.^t{y)- In this notation, the 
objective function is redefined as 

f{y) = '^Pufu{y), (IV.6) 

u=l 


where, for any t = 1,2, ...,T, fu{y) = Ececexp(-ac,t|/c,t)-s“t(|/), and elements 
of the gradient V fu{y) are 


dfuiy) 

dyc,t 




(IV.7) 


The OA algorithm with U cuts per iteration is described as Algorithm 9. 
The algorithm solves the following master problem MP2(/c) in the k-ih. iteration. 
We denote the corresponding optimal value and optimal solution for and y^^\ 
respectively. 

Formulation of Master problem : MP2(/c) 


min ^ = Y!i=i PuZu 
s.t. 

Zu > fuiy^"'^) + V/„(|/(*))'(?/ - !/(*)) M = 1, 2,..., f/, f = 0,1,..., /c - 1 
5I/c'e7^(c) — X]c'eJP(c) ^c,c',t = 1,..., T 

yc,t X/c'eJ^(c) ^C,c',t ^ ^ 

Xc,c',u yc,t- integer 


Algorithm 9 (OA Algorithm with multi-cut). 

Step 0. Set a relative optimality tolerance 5 > 0. Choose an initial feasible 
solution y^^\ 

Step 1. Calculate fu{y^^^) and V/„(|/*'°^), m = 1, 2, ..., U (see (IV.6) and (IV.7)). 
Set lower bound g = 0, upper bound q = Y.u=iPufuiy^^^), and k = 1. 
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Step 2. If the gap {q—q)/q < S, stop. Else, solve the master problem MP2(A;), 
and obtain its optimal value and optimal solution If zi^) > q, then 
q = z^’^\ 

Step 3. Calculate and Vu = 1,2,..., U (see (IV.6) and (IV.7)). 

Step 4. If /(|/(^)) = Y.u=iPufu{y^'''^) < q, then q = f{y^^^). Replace k by 
k + 1, and go to Step 2. 

Similar to the master problem MP1(A;) in Algorithm 8, the relative optimality 
tolerance of the master problem MP2(/c) is denoted by > 0. Again we note that, 
if > 0, z^^^ must be replaced by (1 — to obtain a valid lower bound in 

Step 2 of Algorithm 9. Here, z^^}) is the value turn by the solver when using relative 
optimality tolerance > 0 when solving MP2(A;). 

We examine Algorithm 9 using the same problem instance as in Table 17. In 
these instances, conditioning on the target movement in the initial two time steps 
{t=l and 2) results in U=5 scenarios. Conditioning on the initial three time steps 
results in 7/ = 25 scenarios. We refer to these two versions of Algorithm 9 as “multi- 
t2-cut” and “multi-t3-cut,” respectively. As above, the relative optimality tolerance 

when solving the master problem MP2(/c) is set to 0 or 0.03. This results in a 
total of four versions of Algorithm 9: 

• version 1: multi-t2-cut, and MP2(A;) solved with = 0 

• version 2: multi-t2-cut, and MP2(A;) solved with = 0.03 

• version 3: multi-t3-cut, and MP2(A;) solved with = 0 

• version 4: multi-t3-cut, and MP2(/c) solved with = 0.03 

We note that the initial feasible solution is still set to be the trivial search 
plan (i.e., all searchers keep searching the initial cell for the whole duration). Table 
18 shows the numerical results for versions 1 and 2 of Algorithm 9. Versions 1 and 
2 spend 96.3% and 80.1% of the run time to solve the master problems (in Step 
2), respectively. The remaining time is used to generate cuts and build models. 


86 



During one hour, the version 1 permits only 20 iterations while version 2 executes 
68 iterations. Since the multi-cut requires much more overhead time to generate 
(multiple) cuts than Algorithm 8 (single-cut) in each iteration, the number of possible 
iterations is reduced compared to that of Algorithm 8 (i.e., 25 iterations in version 1 
of Algorithm 8 are reduced to 20 in version 1 of Algorithm 9; 96 iterations in version 
2 of Algorithm 8 are reduced to 68 in version 2 of Algorithm 9). However, for both 
version 1 and version 2 of Algorithm 9, the total nnmber of accnmnlated cuts becomes 
larger than the corresponding versions using a single-cnt approach. As a result, the 
lower bound of the multi-t2-cut (versions 1 and 2) improves qnicker than that of the 
single-cnt (refer to Table 17). For the npper bonnd, there is no signihcant difference 
between the mnlti-t2-cut and the single-cnt approaches. After one honr, the multi- 
t2-cnt versions provides better relative gaps than the single-cnt versions (i.e., 30.4% 
(version 1 of Algorithm 8) compared to 28.3% (version 1 of Algorithm 9), and 20.5% 
(version 2 of Algorithm 8) compared to 19.6% (version 2 of Algorithm 9)). We also 
note that version 2 = 0.03) performs better than the version 1 (e^^^ = 0)) as seen 

from Table 17. 


Time 

Algo. 9 version 1: e 

(fc) = 0 

Algo. 9 version 2: e^) 

= 0.03 

(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.436285 

0.588137 

0.348 

0.470301 

0.569675 

0.211 

20 

0.454153 

0.588137 

0.295 

0.474848 

0.569675 

0.200 

30 

0.457303 

0.588137 

0.286 

0.475066 

0.569675 

0.199 

40 

0.458369 

0.588137 

0.283 

0.475835 

0.569675 

0.197 

50 

0.458369 

0.588137 

0.283 

0.476396 

0.569675 

0.196 

60 

0.458369 

0.588137 

0.283 

0.476396 

0.569675 

0.196 


Table 18. Nnmerical resnlts for Algorithm 9 using mnlti-t2-cnt (versions 1 and 2) 
for every 10 minntes of calculations. The best known non-detection probability is 
0.563472 same as in Table 17. 


Table 19 presents the nnmerical results for Algorithm 9 using mnlti-t3-cnt. In 
each iteration, the mnlti-t3-cut reqnires mnch more time to generate 25 cuts. Thus 
the versions 3 and 4 spend only 40.1% and 4.7% of the rnn time on solving the 
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master problems (in Step 2) and the remaining time generating cuts. During one 
hour tests, versions 3 and 4 permit only 12 and 19 iterations, respectively. Typically, 
the multi-t3-cut versions provide a signihcant improvement in the lower bound in each 
iteration. However, the long run time per iteration prevents them from carrying out 
“enough” iterations. Thus it appears that the multi-t3-cut versions are inferior to the 
single-cut and multi-t2-cut versions (see Tables 17 and 18). We note however that 
our implementation of Algorithm 9 in GAMS could be improved by implementing the 
time-consuming cut generation in C-|--|- or some other faster programming language. 
If the cut generation time could be reduces, the multi-t3-cut versions may prove 
competitive. However, in this dissertation, we do not consider this programming 
enhancement. 


Time 

Algo. 9 version 3: e 

(fc) = 0 

Algo. 9 version 4: 

= 0.03 

(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.255538 

0.711572 

1.785 

0.247878 

0.716509 

1.891 

20 

0.359439 

0.624907 

0.739 

0.343276 

0.662075 

0.929 

30 

0.399356 

0.604164 

0.513 

0.384563 

0.608274 

0.582 

40 

0.424481 

0.604164 

0.423 

0.419195 

0.605206 

0.444 

50 

0.424481 

0.604164 

0.423 

0.435430 

0.589771 

0.354 

60 

0.429744 

0.594532 

0.383 

0.450207 

0.581592 

0.292 


Table 19. Numerical results for Algorithm 9 using multi-t3-cut (versions 3 and 4) 
for every 10 minutes of calculations. The best known non-detection probability is 
0.563472 same as in Table 17. 


Based on the test results in Tables 17, 18 and 19, we conclude that Algorithm 
8 , version 2 (single-cut with e^^^=0.03), and Algorithm 9, version 4 (multi-t2-cut with 
g(*:)=0.03), are the most promising. Hence, we adopt these versions as baselines for 
further comparisons. When the meaning is clear from the context, we refer to those 
versions simply as single-cut and multi-cut, respectively. 

(b) Strong-cut 

In Algorithms 8 and 9, at each iteration, cuts are generated at the optimal so¬ 
lutions y obtained from the master problems. The cuts are the tangent hyper-planes 




y y+i 


Figure 11. Strong cut and tangent hyper-plane on an exponential objective function 
in mixed-integer program. 

and utilize that f{y) > f{y) -|- fiyYiy — y) for all y. (For simplicity, we now argue 
using single-cut, but the same arguments can be build for multi-cut versions). We 
now use the fact we are only concern with integer values of y to build a stronger cut. 
Figure 11 illustrates the basic idea of the new cut. Since the MHSP is an integer 
program, we utilize hnite differences of the objective function f{y) by considering the 
perturbation from yc,i to yc,t + 1 while keep all other variables hxed. Theorem IV. 1 
presents this new cut. Let /\c,t ^ be a vector in which the (c, t) element is one 
and the other elements are all zero, where Z is the set of integers. 
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Theorem IV.l For any y, y e , f{y) > f{y) + Ec,t[f(y+^c,t)-f{y)]{yc,t-yc,t)- 


Proof. f{y) = Ea.Pa.exp(-Ec,t<tl/c,t) = Ea;exp(-Ec,i <tl/c,t+logPa.). Let a‘^ be 
a CT-dimensional vector defined by a‘^ = and = (—logp^). Then, > 0 and 
> 0. Hence, f{y) = fui{y) where = exp(—a‘^1/— 6*^). In this notation, the 
result is equivalent to f^{y) + Ec,J/a;(y + ^c,t) - fUy)]iyc,t - yc,t) < fM for aH a;. 
Consequently, fujiy)[l + Ec,i(exp(-ac,t) “ - yc,t) - exp{-a‘^{y - y))] < 0. Let 

= exp(—and Af be the set of the cell-time pair (c, t) such that yc,t — yc,t 7^ 0 
and > 0 (i.e., cell c is on path u at time period t). Then we only need to show 
that EieAr(l - Pi){yi - m) + > 1, where i = (c,t) G A/”, 0 < A < 1. 

Let ki = yi - jji, then we can define ip{P) = EieAr(l - where 

(3 = (A) is a |A/'|-dimensional vector. Let '^(/3) = HieArA^L Then VipiP) = 
(log(A), •••, fog(AA^|))A(/^)) and V^'^(/5) = where H is a lA/”] by lA/”! matrix 

in which the first low is (log(A),•••, fog(AA^l)) and the other elements are all zero. 
Since the matrix A'A is positive semi-definite and 'ip{j3) > 0, '^(/3) is convex on the 
set with A > 0 for all i G Af. Thus v^(/3) is also convex on the set with A > 0 for 
all i G A/”. When A=1 for Vi G A/”, Vv?(/d) = (A(-l + •••, ^|Ar|(-l + /3|^|)) = 0. 

Thus, the minimum value of 9?(/5) is 1 when A=1 for Vi G Af. Consequently, we 
obtain that 9?(/5) > 1 which proves the result. □ 

We refer to this type of new cut as ’strong-cut’. The calculation of finite 
differences [f{jj + Z.c,t) — f(y)] is as easy as that of the gradient (IV.5). Specifically, 

lf(y + f^c,t) - f(y)] 

= rc',t(y)[exp(-ac',tyc',t - Oic,t) - exp{-ac',tyc',t)\sc',t{y) 

c'ec 

= rc,t(|/)[exp(-Q;c,t(l/c,t + 1)) - exp{-ac,tyc,t)]sc,t{y)- (IV.8) 

For the multi-cut version, for each scenario u = the finite differences 

[fu{y + — fu{y)] are calculated similarly: 
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(IV.9) 


[fu{y + Ac,t) - fu{y)] 

= H ^^c',t(l/)[exp(-ttc',tl/c',i - ac,t) - exp{-ac',tyc',t)K',t{y) 
c'ec 

= rlt{y)[exp{-ac,t{yc,t + 1 )) - exp{-ac,tyc,t)]slt{y)- 

Algorithm 10 describes a single-cut OA algorithm using the strong-cut. Algo¬ 
rithm 10 is a modihcation from Algorithm 8 since at each iteration a strong cut is 
generated instead of a tangent hyper-plane. In the k-th iteration of Algorithm 10, we 
solve a master problem MP3 (A;) dehned below, where the optimal value and optimal 
solution are denoted and y^’^\ respectively. 


Formulation of Master problem : MP3(/c) 


min z 
s.t. 

z > /(i/») + UcAfiy^^ + Ae,t) - /(i/»)](i/c,i - y^l) * = 0,1,..., k - 1 

5Zc'e7^(c) ^c',c,t—1 X/c'eJ^(c) Vt 1, ...,T 

yc,t X/c'eJP(c) ^c,c',t Vt G T 
Xc,c',t, yc,t- integer 


Algorithm 10 (OA Algorithm with single, strong-cut). 

Step 0. Set a relative optimality tolerance 5 > 0. Choose an initial feasible 
solution y^^\ 

Step 1. Calculate f{y^^^) and -|- Ac,t) — f{y^^^) for V(c, t) (see (IV.4) 

and (IV.8)). Set lower bound q = 0, upper bound q = /(|/*'°^), and A; = 1. 

Step 2. If the gap {q—q)/q < S, stop. Else, solve the master problem MP3(A;), 
and obtain its optimal value z^^^ and optimal solution y^^\ If z^^'^ > q, then 
q = z^^\ 

Step 3. Calculate and for y{c,t) (see (IV.4) 

and (IV.8)). 

Step 4. If f{y^’^^) < g, then q = f{y^^^). Replace A; by A; + 1, and go to Step 

2 . 
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Algorithm 11 shows the multi-cut version of Algorithm 10. In the k-ih. it¬ 
eration of Algorithm 11, the master problem MP4(A;) defined below is solved. The 
optimal value and optimal solution of MP4(A;) are denoted and respectively. 
In each iteration of Algorithm 11, U cuts are generated at once. 


Formulation of Master problem : MP4(/c) 


min z = Y!i=i PuZu 
s.t. 

Zu > fu{y^^) + + Ac,i) - fu{y^"'^)]{yc,t - yc}) 

n = 1,2,..., f/, i = 0,1,...,/c — 1 

Ec'e7?.(c) ^c',c,t—1 Ec'eJP(c) Vt 1, ...,T 

yc,t Ec'eJ^(c) ^c,c',t Vf G T 
Xc,c',t, yc,t- integer 


Algorithm 11 (OA Algorithm with multi/strong-cut). 

Step 0. Set a relative optimality tolerance <5 > 0. Choose an initial feasible 
solution y^^\ 

Step 1. Calculate fu{y^'^^) and fu{y^^^ + Ac,t) - fu{y^^^), u = 1,2,...,U, V(c, t) 

(see (IV.6) and (IV.9)). Set lower bound g = 0, upper bound q = En=i Pufu{y^^^), 
and k = 1. 

Step 2. If the gap {q—q)/q < S, stop. Else, solve the master problem MP4(A;), 
and obtain its optimal value z^^^ and optimal solution y^^\ If z(^) > q, then 
q = z^^\ 

Step 3. Calculate /„(|/(^)) and /„(|/(^) -h A^,*) -/„(?/('')), n = 1,2, ...,C, V(c,t) 

(see (IV.6) and (IV.9)). 

Step 4. If f{y^^'>) = Y!i=iPufu{y^'"'^) < q, then q = f{y^^^). Replace k by 
k + 1, and go to Step 2. 

We demonstrate the effect of the strong-cut using the same problem instance 
accompanied with the same assumptions and parameter settings as in the previous 
subsection. The trivial search plan is still used as the initial feasible solution y^'^k 
The relative optimality tolerance of the master problems is assumed constant at 0.03. 
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Table 20 describes the effect of the strong-cut. Compared with the previous 
algorithms (Algorithms 8 and 9), Algorithms 10 and 11 implement 1.4 times more 
iterations: 96 for Algorithm 8 compared to 135 for Algorithm 10, and 68 for Algorithm 
9 compared to 95 for Algorithm 11. Furthermore, the relative gap after one hour is 
drastically reduced from about 20% to 10.5% for both versions. Since the increase 
of the number of iterations and the progress in the gap are essentially identical for 
both versions, we conclude that the strong-cut is quite effective. The performance of 
Algorithms 10 and 11 during one hour of calculations is almost parallel to each other. 


Time 

Algo. 10 : single-cut 

Algo. 

11 : multi-cut 

(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.509594 

0.570631 

0.120 

0.510988 

0.571003 

0.117 

20 

0.511369 

0.570631 

0.116 

0.512539 

0.569664 

0.111 

30 

0.512452 

0.568363 

0.109 

0.513929 

0.569664 

0.108 

40 

0.512819 

0.568363 

0.108 

0.514792 

0.569664 

0.107 

50 

0.512847 

0.568363 

0.108 

0.515407 

0.569664 

0.105 

60 

0.514277 

0.568363 

0.105 

0.515407 

0.569664 

0.105 


Table 20. Numerical results for Algorithms 10 and 11 using strong-cut for every 10 
minutes of calculations. The best known non-detection probability is 0.563472 same 
as in Table 17. 


(c) Choice of Initial Solution 

So far, as a choice of initial feasible solution, we have used the trivial search 
plan (i.e., all searchers keep loitering at the initial cell). The choice of the initial solu¬ 
tion certainly influences the OA algorithms, especially during the initial calculations. 
In this section, we aim to develop a “good” initial solution. In advance of this, we 
hrst introduce a theoretical property of the MHSP, which facilitates the calculation 
of a good initial solution. 

We now assume that the detection rate is time and space independent, i.e., 
a = ac,t for all c and t. Under this assumption, we show that the MHSP is equivalent 
to a large-scale linear mixed-integer program. Let be the total 

effective search effort given search plan y when the target takes path u & Q. Since 
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the detection rate is assumed identical, is a multiple of a and thus = ma for 
some m = 0,1,JT. This is a similar idea to [52], The maximum possible total 
effective search effort is the number of searchers (J) times the time horizon (T), which 
occurs when all searchers take the path oj. In this notation, we claim that the MHSP 
is formulated as the following large-scale linear mixed-integer program. 

Formulation in linear mixed-integer program: LI 


Indices 

c, c' cells (c, c' G C = { 1 ,..., C}) 
t time step (f G T = {0,1,..., T}) 

j searchers {j E J = {1,..., J}) 

OJ target path (a; G ff) 

Parameters 

a detection rate for a single searcher in cell c and time period t 

ac,t if cell c is on target path u at time period t, zero otherwise. 
i/cfl number of searchers in cell c in time period 0 
probability that the target takes path oo 

Variables 

Xc^c't number of searchers that is redistributed from cell c in time 
period t to cell c' in time period t + 1 
yc,t number of searchers in cell c in time period t 

Formulation 

min 

s.t. 

Zu) > -f m — me““] -|- (1/Q;)e“'^"(e“" — l)Wui 

Woo G n, m = 0,1..., JT 
5Zc'e7^(c) 1 ^c'£j^{c) ^C,c',t 'Jt 1, ...,T 

Ucj X/c'eJ^(c) 'Jt G T 

Xc,c',t, VcT integer 

Theorem IV. 2 If the detection rate is identical (i.e., a = ac,t for all c and t), the 
convex nonlinear mixed-integer MHSP and the linear mixed-integer program LI have 
identical global minimum solutions. 

Proof. Let UiWu;) = where = Ec,t and f(g) = Eu;eQPcufuj(W^)- 

only takes a hnite number of discrete value ffE = ma, m = 0,1, ..., JT. Thus 
fojiWoj) can be described as a piece-wise linear function (see Figure 12). The linear 
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Figure 12. Piece-wise linearlization of the exponential function fcj{W^). 

approximation between ma and (m -|- 1 )q; is obtained as f(W^) = + m — 

me~°‘] + (1/Q;)e“™“(e““ — 1)W^ by solving the linear equation {bo + 6i[mQ;] = 
and bo + bi[{m + l)a] = for parameters bo and bi. □ 

The linearization of MHSP under the assumption of identical detection rates 
explicitly depends on all possible target paths. Hence, the program LI may become 
extremely large. Thus we cannot directly apply this formulation to solve most in¬ 
stances of MHSP. However, this linearization provides useful insight to obtain a good 
initial solution. 

The linearlization of MHSP considers the total effective search effort IT^ on 
each target path cu G H individually. The alternative approach is to consider the 
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total effective search effort aggregated over all target paths. This approach allow us 
to utilize the Markovian property of the target movement, but it leads to a relaxation 
as we now described. Recall that the target moves between cells according to a 
known transition matrix T. Let p{-,t) = t > 1 be the undetected target 

distribution in time period t, (p(-, 1) is known). Thus the aggregated total effective 
search effort is described as W = J2c,tP{c,t)ayc,y In this notation, the aggregation 
of the large-scale linear mixed-integer program LI takes the form of the moderately 
sized linear mixed-integer program L2: 

Formulation in linear mixed-integer program (aggregation): L2 


Indices 

c, c' cells (c, c' G C = (1,..., C}) 
t time step (t G T = (0,1,..., T}) 

j searchers {j E J = {1,..., J}) 

Parameters 

a detection rate for a single searcher in any cell and time 
i/cfl number of searchers in cell c in time period 0 
p{-,t) undetected target distribution at time period t 

Variables 

Xc^c't number of searchers that is redistributed from cell c in time 
period t to cell d in time period t + 1 
yc,t number of searchers in cell c in time period t 

Formulation 

min z 

s.t. 

W = Ec,tP(c, t)ayc,y 

z > + m — me““] -1- (l/a)e“'""(e“" — l)W, m = 0,1..., JT 

Ec'&TZ{c) ^c' X)c'eJ^(c) 'dt 1, ...,T 

yc,t Ec'GT{c) 'dt G T 

Xc,c',t, Vcp integer 

We observe that solutions of L2 is not identical to solutions of LI as L2 is a 
relaxation of LI due to the aggregation of all the target paths. However, L2 is easily 
solved and typically provides a good initial solution for the main calculations. We 
will use this aggregated program L2 to obtain an initial solution. 
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In addition, we also consider a simple condition (constraint) to obtain a better 
initial solntion. If the nnmber of searchers is moderate (e.g., 3 searchers) and the size 
of the area is large, it is typically better that the searchers visit cells not previously 
searched. (Of course, the effect of this strategy also depends on the movement of the 
target.) Given a specihc G T (e.g., t'=l for the problem instance with 3 searchers), 
a non-revisit constraint is described as yc,t < 1, Vc G C. Therefore, in the initial 
step (Step 0) of the OA algorithms, we solve the aggregation problem (L2) with this 
non-revisit constraint to obtain an initial solution. This initial problem is referred as 
L3 as follows. 

Initial problem: L3 


Indices 

c, c' cells (c, c' G C = {1,..., G}) 
t time step (t G T = {0,1,..., T}) 

j searchers {j E J = {1,..., J}) 

Parameters 

a detection rate for a single searcher in any cell and time 
Ucfi number of searchers in cell c in time period 0 
p(-,t) undetected target distribution at time period t 

Variables 

Xc,c't number of searchers that is redistributed from cell c in time 
period t to cell d in time period t + I 
Uc^t number of searchers in cell c in time period t 

Formnlation 

min z 

s.t. 

^ = Ec,tP(c, t)ayc,y 

z > + m — me~°‘] + (1/Q;)e“'™"(e“" — 1)1T, m = 0,1..., JT 

5I/c'e7^(c) ^c'— X]c'ej^(c) dt = 1,..., T 
yc,t X)c'eJP(c) dt G T 

J2t>t' yc,t < 1, Vc G C 
yc,t- integer 


In the OA algorithms, the hrst cut(s) is generated at the initial solution ob¬ 
tained from the initial problem L3. In addition, we also add the cut at the solution 
y = 0 (i.e., f{y) > /(O) -I- V/(0)'|/) which can be shown to relate to the mean bound 
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in a branch-and-bound procedure [34, 52], The cut at y = 0 linearlizes the probability 
of non-detection at the point of “no search” effect so it gives an estimate of initial 
effect of increasing the search effort above zero in a cell. We apply the strong-cut, 
not the gradient-based hyper plane, at y = 0 and the initial solution found when 
solving L3. We note that this initial process is also available for the algorithm using 
multi-cut (Algorithm 11) and is referred as the initial problem L4 as follows. 

Initial problem: L4 


Indices 

c, d cells (c, c' G C = { 1 ,..., C}) 
t time step (t G T = {0,1,..., T}) 

j searchers {j E J = {1,..., J}) 

u scenario (m = 1 , 2 ,..., f/) 

Parameters 

a detection rate for a single searcher in any cell and time 
Ucfl number of searchers in cell c in time period 0 
pu probability that the target takes scenario u 

undetected target distribution at time period t when scenario u 

Variables 

Xc^c't number of searchers that is redistributed from cell c in time 
period t to cell d in time period t + 1 
Uc^t number of searchers in cell c in time period t 
Wu total effective search effort when scenario u 

Formnlation 

min YZ=iPuZu 

s.t. 

fb’n = Ec,tP“(c,t)a|/c,y, M=l, 2 ,...,f/ 

Zn > e-™“[l + m - me-“] + (l/a)e-™“(e-“ - l)Wu 

u = 1,2, m = 0, 1..., JT 
Z)c'e 7 ^(c) Xc\c,t-l = J2c'£r{c) Xc,c',t dt = 1,..., T 
yc,t X)c'eJP(c) dt G T 

Et>t’yc,t < 1, Vc G C 
Xc,c',t, yc,t- integer 

Table 21 reports the effect of the choice of an improved initial solution using 
the same problem instance with the same assumptions and parameter settings as 
above. We dehne Algorithm 10.2 (single-cut) and Algorithm 11.2 (multi-cut) when 
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an initial solution is computed from L3 in Algorithm 10 and from L4 in Algorithm 
11, respectively. For both Algorithms 10.2 and 11.2, the relative gap decreases after 
10 minutes from 0.107 to 0.103 for Algorithm 10.2, and from 0.105 to 0.098 for 
Algorithm 11.2. However, after 30 minutes of calculations, the advantage of the 
improved initial solution becomes insignihcant since the numbers of cuts are added 
and the cuts largely influence the solution quality. Algorithms 10.2 and 11.2 perform 
better than the previous cases using trivial initial solution, thus we apply them as 
the baseline for the OA algorithms which include the next set of new cuts. 


Time 

(min) 

Algo. 

LB 

10.2 : single-cut 

UB Gap 

Algo. 

LB 

11.2 : multi-cut 

UB Gap 

10 

0.510520 

0.568875 

0.114 

0.512064 

0.567808 

0.109 

20 

0.512370 

0.568875 

0.110 

0.513271 

0.567808 

0.106 

30 

0.512422 

0.568875 

0.110 

0.514057 

0.567808 

0.105 

40 

0.513546 

0.568875 

0.108 

0.514660 

0.567808 

0.103 

50 

0.513747 

0.568875 

0.107 

0.515650 

0.567808 

0.101 

60 

0.513747 

0.568875 

0.107 

0.515900 

0.567808 

0.101 


Table 21. Numerical results for Algorithms 10.2 and 11.2 using an improved initial so¬ 
lution for every 10 minutes of calculations. The best known non-detection probability 
is 0.563472 same as in Table 17. 


(d) Relative-cut 

In Table 21, the upper bound (i.e., the best feasible value) does not change 
after ten minutes. Thus it is possible that the upper bound could have been optimal. 
(Of course, it is not in this case as we know a better value from Chapter III). Let y 
be the best solution so far providing the upper bound q = f{y). Since the objective 
function is convex, y is guaranteed to be optimal if f{y) is still the upper bound after 
cuts are added at all neighbor integer points y oi y (the Euclidian distance between 
y and y is one). Thus, we explore the strategy of generating cuts around the best 
solution. 

Figure 13 illustrates the rule that determines a neighboring integer point y 
where we generate a cut. In earlier the OA algorithms, each iteration obtains a 
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Figure 13. Illustration of cut generation around the best solution. 

solution y and generates one or more cuts at y. Now, we also generate an additional 
cut at y. Thus, in the algorithms using the single-cut, two cuts are generated in each 
iteration. The upper picture in Figure 13 shows a case where the best solution so far 
y is still better than the current solution y (i.e., f{y) < f{y)). Then, we determin 
an integer point y which is a neighbor of y and is the closest point to y measured in 
the Euclidean distance. The lower picture illustrates the opposite case. In this case, 
after generating the cuts at y and y, the best solution y is updated as y. We denote 
the cut at the solution y as a “relative-cut. ” Moreover, for the relative-cut, we also 
apply the strong-cut instead of the gradient-based hyper-plane cut. 

The OA algorithms using the relative-cut is described as Algorithm 12 and 
Algorithm 13. Algorithm 12 is the single-cut version, and Algorithm 13 is the multi- 
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cut version. Algorithms 12 and 13 also use the initial problem L3 and L4 to obtain a 
reasonable initial solution, respectively. Thus Algorithms 12 and 13 are extension of 
Algorithms 10.2 and 11.2 to the ones with relative-cut, respectively. 

Algorithm 12 (OA Algorithm with single/strong/relative-cut). 

Step 0. Set a relative optimality tolerance 5 > 0 and = 0. Solve the 
initial problem L3 {t' = 1 for the problem instance with 3 searchers) and 
set the solution as 

Step 1. Calculate /(|/(°^), f{y^^'> + Ac,t) - f{y^^'>) and f{y^^'^ + /\c,t) - 

f{y^^'^) for V(c,t) (see (IV.4) and (IV.8)). Set lower bound g = 0, upper 
bound q = y = y^^\ and k = 2. 

Step 2. If the gap {q — q)/q < S, stop. Else, solve the master problem MP3(/c) 
(described above Algorithm 10), and obtain its optimal value and 
optimal solution y = y^^\ If > g, then q = z^^\ 

Step 3. Calculate and + /2zc,t) — for V(c, t) (see (IV.4) 

and (IV.8)). 

Step 4. Obtain a neighbor point y around the best solution by comparing 
f{y) with q = f{y). Replace k hj k + 1. Set = y. 

Step 5. Calculate and + /2zc,t) — for V(c, t). 

Step 6. If f{y) < q, then q = f{y) and y = y. Replace k hy k + 1, and go to 

Step 2. 

Algorithm 13 (OA Algorithm with multi/strong/relative-cut). 

Step 0. Set a relative optimality tolerance 5 > 0 and = 0. Solve the 
initial problem L4 [f = 1 for the problem instance with 3 searchers) and 
set the solution as y^^\ 

Step 1. Calculate /«(|/^°^), fuiy^^'^), Uiy^^^ + ^c,t) - /«(l/^°^) and fuiy^^'^ + 

^c,t) — fuiy^^^), for n = 1,2,...,!/, V(c, t) (see (IV.6) and (IV.9)). Set lower 

bound g = 0, upper bound g = Y.u=iPufu{y^^^), y = y‘^^\ and k = 2. 

Step 2. If the gap (g — g)/g < 5, stop. Else, solve the master problem MP4(/c) 
(described above Algorithm 11), and obtain its optimal value and 
optimal solution y = y^^\ If z^^'^ > g, then g = z'^^k 
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step 3. Calculateand+ A^,*)u = 1,2, ...,A, V(c,t) 

(see (IV.6) and (IV.9)). Calculate f{y) = T,^=iPufu{y) 

Step 4. Obtain a neighbor point y around the best solution by comparing 
f{y) with q = f{y). Replace k hj k + 1. Set = y. 

Step 5. Calculateand+ n = 1,2, ...,C, V(c,t). 

Step 6. If fijj) < q, then q = f{y) and y = y. Replace k hy k + 1, and go to 
Step 2. 

We test the effect of the relative-cut using the same problem instance with 
the same assumptions and parameter settings as above. Table 22 describes the com¬ 
putational results of Algorithms 12 and 13. During one hour of calculations, we do 
not identify signihcant effect of the relative-cut by comparing Table 21 and Table 22. 
However, the relative-cut strategy seems not to be detrimental and it allows for a 
faster accumulation of cuts which may be benehcial so we implement the relative-cut 
in our OA algorithms. 


Time 

Algo. 

12: single-cut 

Algo. 

13: multi-cut 

(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.510523 

0.571749 

0.120 

0.512066 

0.577300 

0.127 

20 

0.512079 

0.569864 

0.113 

0.514172 

0.570003 

0.109 

30 

0.513185 

0.569864 

0.110 

0.514721 

0.566784 

0.101 

40 

0.514212 

0.569864 

0.108 

0.514806 

0.566784 

0.101 

50 

0.514538 

0.569864 

0.108 

0.515639 

0.566784 

0.099 

60 

0.514538 

0.569864 

0.108 

0.515690 

0.566784 

0.099 


Table 22. Numerical results for Algorithms 12 and 13 using relative-cut for every 10 
minutes of calculations. The best known non-detection probability is 0.563472 same 
as in Table 17. 


(e) Symmetric-cut 

The hnal cut we present is quite specihc and not applicable for general cases. 
However it is somewhat effective in specihc situation. Suppose that the shape of the 
area of interest (AOI) is symmetric and the searchers and the targets are initially in 
the cells which are located on a line that divides the AOI symmetrically (see Figure 
14). Suppose also that the detection rate is identical for all cells and time periods. 
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And finally, suppose that the target moves to an adjacent cells equally likely. We 
refer to this specihc situation as “symmetric.” In the symmetric situation, for a 
search plan V, there exist a symmetric (or mirror) plan V' with identical probability 
of nondetection, i.e., q{V) = q{V'), see Figure 14. 

Thus, the basic idea of the symmetric-cut is that, in each iteration, cuts are 
generated not only at the current solution (plan) V but also at its symmetric solution 
(plan) v' since the symmetric solution can be evaluated from the current solution. 
Algorithm 14 describe the OA algorithm using the symmetric-cut, which is an exten¬ 
sion of Algorithm 12. In Algorithm 12, at each iteration, two cuts are generated at 
the current solution y and a neighbor point y around the best solution (by comparing 
fijj) with fiy)). In Algorithm 14, two cuts are additionally generated at the sym¬ 
metric solution y' and a neighbor point y' around the best solution (by comparing 
f{y') with f{y))- Thus, in Algorithm 14, four cuts are added at each iteration. 

The problem instance examined in numerical test so far is symmetric. Thus 
we attempt to examine the effect of the symmetric cut in this instance. Since the 
multi-cut utilizes the conditioning of the target movement, it violates the symmetry 
property. Thus the symmetric cut is available in the algorithm using the single-cut 
only. 

Algorithm 14 (OA Algorithm with single/strong/relative/symmetric-cut). 

Step 0. Set a relative optimality tolerance 5 > 0 and y^^^ = 0. Solve the initial 
problem L3 (described below Algorithm 12, and t' = 1 for the problem 
instance with 3 searchers) and set the solution as y^^\ 

Step 1. Calculate f{y^°^), f{y^^'> + A^,*) - f{y^^'>) and f{y^^'^ + /\c,t) - 

f{y^^^) for V(c,t) (see (IV.4) and (IV.8)). Set lower bound q = 0, upper 
bound q = f{y^^^), y = y^^\ and k = 2. 

Step 2. If the gap {q — q)/q < S, stop. Else, solve the master problem MP3(/c) 
(described above Algorithm 10), and obtain its optimal value and 
optimal solution y = y^^\ If > g, then q = z^^\ 
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Search plan P 
Symmetric plan P 





Figure 14. Symmetric environment where the symmetric-cut is available. 

Step 3. Calculate and -|- Ac,t) — for V(c, t) (see (IV.4) 

and (IV.8)). 

Step 4. Obtain a neighbor point y around the best solution by comparing 
f{y) with q = f{y). Replace k hj k + 1. Set = y. 

Step 5. Calculate and -|- Ac,t) — f{y^^^) for V(c, t)- 

Step 6. For the solution y, obtain the symmetric solution y'. Replace k by 

k + 1. Set = y'. 

Step 7. Calculate and -|- Ac,t) — f{y^^^) for V(c, t)- 

Step 8. Obtain a neighbor point y' around the best solution by comparing 
f{y') with q = f{y). Replace khy k + 1. Set y^'^^ = y'. 

Step 9. Calculate and -|- /\c,t) — f{y^^^) for V(c, t)- 

Step 10. If f{y) < q, then q = f{jj) and y = y. Replace khy k + 1 and go to 
Step 2. 
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Table 23 reports the effect of the symmetric-cut as applied to the same problem 
instance as in Table 22. Comparing with Algorithm 12 in Table 22, the relative gap 
after one honr of calculations is reduced from 0.097 to 0.092. Both the lower and 
upper bound are improved. Thus, among the OA algorithms using the single-cnt. 
Algorithm 14 is the best for this problem instance. 


Time 

Algo. 14: 

Symmetric-cut 

(min) 

LB 

UB 

Gap 

10 

0.512580 

0.574371 

0.121 

20 

0.514023 

0.572793 

0.114 

30 

0.515390 

0.568628 

0.103 

40 

0.515390 

0.568628 

0.103 

50 

0.515812 

0.568628 

0.102 

60 

0.516447 

0.568628 

0.101 


Table 23. Nnmerical results for Algorithms 14 using symmetric-cnt for every 10 
minntes of calculations. The best known non-detection probability is 0.563472 same 
as in Table 17. 


We examined several OA algorithms in which the new cnts (multi-, strong-, 
relative- and symmetric-cnts) are cnmnlatively applied. Specihcally, Algorithm 13 
(multi-cnt) and Algorithm 14 (single-cnt) perform best among the algorithms we 
examined. For the problem instance (15 by 15 cells, time horizon 18 and 3 searchers) 
with specific assnmptions and parameters settings (e.g., identical detection rate and 
symmetric situation, etc.). Algorithms 13 and 14 provide an approximate solution 
with relative optimality gap about 10% after one hour of calculations. In Chapter 
III, for the same problem instance (see Tables 15 and 16), the Cross-Entropy (CE) 
henristic provides a better solntion 0.563472 (=1-0.436528) after 674.85 seconds. Thus 
the CE heuristic performs better for this one-hour computational test, bnt, of course 
the CE heuristic does not provide solution qnality guarantees. Thus a hybrid method 
using the CE henristic and the OA algorithm is expected to be more effective: the 
CE heuristic is initially applied to obtain a good initial feasible solntion followed by 
iterations of the OA algorithm. This dissertation does not pursne this hybrid scheme. 
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3. Numerical Study of Three Searchers 

In the previous subsection, we developed specific OA algorithms (i.e., Algo¬ 
rithm 13 and Algorithm 14 among others) to solve the MHSP. For the moderately 
sized problem instance (15 by 15 cells, time horizon 18 and 3 searchers) examined, 
the OA algorithms provided a solution about 10% of the optimal solution in one hour. 
We now examine the run time of Algorithm 13 on a more extensive set of problem 
instances of somewhat smaller size. Results for Algorithm 14 is not presented, but 
they are almost identical. 

As problem instances, we consider the following five cases with three searchers 
in which we apply the same assumptions and parameter settings as the problem 
instance in the previous subsection, except for changing the area size and the time 
horizon. 


• case 1: 

• case 2: 

• case 3: 

• case 4: 

• case 5: 


15 by 15 cells and time horizon 16 
13 by 13 cells and time horizon 14 
11 by 11 cells and time horizon 12 
9 by 9 cells and time horizon 10 
7 by 7 cells and time horizon 8 


For each case, we run Algorithm 13 for two hours. Since the problem instances 
are small, the master problems MP4(/c) in Algorithm 13 are expected to be solved 
quickly. Thus, we allow a smaller relative optimality tolerance for the solver of 
MP4(A;). Specifically, we find empirically that the following simple rule works well; 
set 6*^^^= min{0.03,(g — 5')/(3g)}. In short, at each iteration, we use the optimality 
tolerance 0.03, but apply one third of the latest gap (g — q)/q after the gap has 
dropped below 9%. 

Table 24 reports the gaps at every ten minute during two-hours numerical tests 
using Algorithm 13 on these five cases. Essentially all cases obtain solution within 
5% of optimality within 20 minutes. The less the size of the problem instance is. 
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the smaller gaps are achieved. Case 5 results in an essentially optimal solution after 
two hours. After two hours, the lower and upper bounds in that case are 0.454141 
and 0.455054, respectively. Thus Algorithm 13 is quite effective for these problem 
instances and provides near-optimal solutions in reasonable computing times. 


Gap 

easel case2 case3 easel ease5 


Time 

(min) 

To 

20 

30 

40 

50 

60 

70 

80 

90 

100 

no 

120 


0.056 

0.051 

0.049 

0.048 

0.048 

0.047 

0.046 

0.044 

0.041 

0.041 

0.041 

0.041 


0.050 

0.047 

0.044 

0.041 

0.040 

0.038 

0.038 

0.037 

0.037 

0.036 

0.036 

0.036 


0.039 

0.034 

0.031 

0.029 

0.028 

0.027 

0.026 

0.026 

0.026 

0.025 

0.025 

0.024 


0.025 

0.020 

0.018 

0.016 

0.015 

0.014 

0.013 

0.012 

0.012 

0.011 

0.010 

0.009 


0.022 

0.016 

0.013 

0.011 

0.009 

0.007 

0.006 

0.005 

0.004 

0.004 

0.003 

0.002 


Table 24. Relative optimality gaps for Algorithm 13 applied three searcher problem 
instances during two hours of calculations. 


4. Application with Large Number of Searchers 

This subsection considers the MHSP with more searchers, especially 5, 10, 15, 
and 30 searchers. Conceptually, the OA algorithms can treat the HMSP with a large 
number of searchers more easily than the Cross-Entropy (CE) heuristics since they 
do not require the generation of specific network structures. Moreover, discussed and 
demonstrated below, the OA algorithms actually may beneht from more searchers 
as the continuous relaxation of the master problems may becomes stronger and thus 
reduce the master problem solution times. In this subsection, we examine the effect 
on the OA algorithms with respect to the number of searchers. For the computational 
tests, we use the same problem instance (15 by 15 cells and time horizon 18) as in 
the previous subsection 2, except we change the number of searchers. The other 
assumptions and the parameter settings are identical except that: (i) The identical 
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detection rate a is adjusted according to the number of searchers. For example, for the 
case with 30 searchers, detection rate O.la is used to make the search effect equivalent 
to the 3 searcher case (i.e., 3 q; = 30 ■ 0.1a). (ii) The optimality tolerance for the 
master problem solver is adjusted adaptively by setting min{0.03,(g — q)/{^q)}- 
For the cases with a large number of searchers, the master problem (mixed-integer 
program: MIP) becomes almost equivalent to the continuous relaxation of the MIP. 
Thus, the master problem with small tolerance (e.g., 0.01) can be solved quickly thus 
we do not need a large optimality tolerance (0.03). On the other hand, the MIP with 
much smaller optimality tolerance (e.g., 0.001) is quite time-consuming to be solved. 
(3) In the initial problem to hnd a good initial solution, the non-revisit constraint is 
ignored since that constraint may be detrimental in situations with a large numbers 
of searchers in a small and moderately sized areas. 

We examine the following four cases are examined on the problem instance 
with 15 by 15 cells and time horizon 18: 

• case 1: 5 searchers 

• case 2: 10 searchers 

• case 3: 15 searchers 

• case 4: 30 searchers 

Since the number of searchers is fairly large, a continuous relaxation of the 
MHSP becomes almost equivalent to the original integer MHSP problem. A con¬ 
tinuous relaxation means that the integrality restriction in MHSP is relaxed, which 
implies that the searchers can be “partitioned” arbitrarily. We denote the instances 
corresponding to cases 1-4 with the continuous relaxation for cases lr-4r. 

For the computational tests, we use Algorithm 13 (multi-cut version) but 
Algorithm 14 generates almost identical results. We refer to Algorithm 13 without 
integer restrictions in the master problem for Algorithm 15. Clearly, Algorithm 15 
may obtain non-integer solutions and thus a search plan that is not implementable. 
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One way of obtaining an (integer) search plan from the non-integer solution is by 
applying some rounding heuristic. However, we have not examined that approach to 
obtain integer solutions from the continuous relaxation. Algorithm 15 does provide 
a lower bound on the optimal value of MHSP as it solves a relaxation. We note that 
gradient-descent methods may also solve the relaxed MHSP. However, we have not 
pursued that avenue. 

Tables 25 to 28 report the results of one-hour computational test for the case 
with 5, 10, 15 and 30 searchers, respectively. For the cases with more than 5 searchers, 
after one hour. Algorithm 13 provides quite good solutions with gaps at about 2%. 
Even for 5 searchers, the gap is moderate about 5%. Algorithm 15 (continuous relax¬ 
ation) gives non-detection probability (NDP) values close to those from Algorithm 
13. Especially for the cases with large number of searchers (15, 30), the obtained 
NDP values from Algorithms 13 and 15 are almost same. 

Comparing with the lower bounds of Algorithms 13 and 15, we hnd that 
Algorithm 15 provides larger lower bounds than Algorithms 13. Of course, the upper 
bound from Algorithm 15 is not valid for MHSP. Algorithm 15 solves each iteration 
more quickly than Algorithms 13 since it solves continuous relaxation problems. Thus, 
Algorithm 15 generates many more cuts and brings up the lower bound quickly. Thus, 
if the allowed calculation time is small, a hybrid method using Algorithms 13 and 
15 may be better. For example. Algorithm 15 (continuous relaxation) is initially 
implemented for some period to push up the lower bound quickly by generating 
many cuts, and later Algorithm 13 is implemented to provide a good (integer) search 
plan. We note the Algorithm 13 with many cuts takes much run time. Thus, this 
hybrid approach may be extremely inefficient when Algorithm 15 has generated too 
many “weak” cuts before Algorithm 13 starts. We examined this hybrid approach for 
several cases, but found it generally to be an inferior approach for this reason. 
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Time 

Algo. 13: case 

1 

Algo. 

15: case 

Ir 

(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.505386 

0.554756 

0.098 

0.506761 

0.547174 

0.080 

20 

0.514396 

0.546065 

0.062 

0.518175 

0.538696 

0.040 

30 

0.516137 

0.546065 

0.058 

0.521317 

0.535791 

0.028 

40 

0.518325 

0.546065 

0.054 

0.523267 

0.534307 

0.021 

50 

0.518325 

0.546065 

0.054 

0.524436 

0.533078 

0.016 

60 

0.518810 

0.546065 

0.053 

0.525279 

0.532586 

0.014 


Table 25. Numerical results for Algorithms 13 and 15 for the case with 5 searchers 
for every 10 minutes of calculations. 


Time 

(min) 

Algo. 13: case2 
LB UB 

Gap 

Algo. 

LB 

15: case2r 
UB 

Gap 

10 

0.501388 

0.549834 

0.097 

0.506429 

0.546805 

0.080 

20 

0.511533 

0.541721 

0.059 

0.518083 

0.538787 

0.040 

30 

0.518910 

0.535177 

0.031 

0.521340 

0.535598 

0.027 

40 

0.520400 

0.535177 

0.028 

0.523113 

0.534203 

0.021 

50 

0.521077 

0.535081 

0.027 

0.524351 

0.533240 

0.017 

60 

0.522048 

0.533785 

0.022 

0.525165 

0.532440 

0.014 


Table 26. Numerical results for Algorithms 13 and 15 for the case with 10 searchers 
for every 10 minutes of calculations. 


Time 

(min) 

Algo. 13: case3 
LB UB 

Gap 

Algo. 

LB 

15: case3r 
UB 

Gap 

10 

0.495162 

0.553566 

0.118 

0.507577 

0.546699 

0.077 

20 

0.511244 

0.540650 

0.058 

0.517370 

0.538297 

0.040 

30 

0.517330 

0.537061 

0.038 

0.520760 

0.536154 

0.030 

40 

0.519869 

0.535204 

0.029 

0.522931 

0.534244 

0.022 

50 

0.521624 

0.534385 

0.024 

0.524083 

0.533416 

0.018 

60 

0.523123 

0.533276 

0.019 

0.524873 

0.532910 

0.015 


Table 27. Numerical results for Algorithms 13 and 15 for the case with 15 searchers 
for every 10 minutes of calculations. 
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Time I Algo. 13: case4 I Algo. 15: case4r 


(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.487626 

0.555553 

0.139 

0.508092 

0.542999 

0.069 

20 

0.505637 

0.542030 

0.072 

0.518134 

0.538112 

0.039 

30 

0.513595 

0.538593 

0.049 

0.521319 

0.535204 

0.027 

40 

0.516808 

0.537237 

0.040 

0.523258 

0.533536 

0.020 

50 

0.520377 

0.534938 

0.028 

0.524341 

0.533536 

0.018 

60 

0.522040 

0.533946 

0.023 

0.524980 

0.532239 

0.014 


Table 28. Numerical results for Algorithms 13 and 15 for the case with 30 searchers 
for every 10 minutes of calculations. 
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V. CONCLUSIONS AND FUTURE 

RESEARCH 

A. CONCLUSIONS 

This dissertation develops models and solution methodologies to solve the 
discrete-time path-optimization problems for single and multiple searchers looking 
for a non-evading moving target in a hnite set of cells. We especially focus on the 
following problems: 

• Single searcher problem (SSP) with additional resource constraints related to 

risk exposure to threats and fuel consumption (RSSP) 

• SSP for multiple searchers (MSP) 

• MSP for multiple homogeneous searchers (MHSP) 

The dissertation starts by formulating RSSP, which generalizes existing models of the 
SSP by considering (i) history-dependent glimpse detection probability, (ii) multiple 
altitudes for the searcher, and (iii) multiple constraints on “consumption” of resources 
such as time, fuel, and risk. 

We develop a specialized branch-and-bound (B/B) algorithm for solving RSSP 
and propose a new bound (Lagrangian directional static bound) on the optimal de¬ 
tection probability using network expansion to account for a portion of the history 
of the current path and using a Lagrangian relaxation to eliminate resource con¬ 
straints. We also derive a series of network reduction procedures that tighten the 
Lagrangian relaxation and reduce the amount of path enumeration required by the 
B/B algorithm. 

In direct comparison with a state-of-the-art algorithm for SSP, the proposed 
bound and network reduction procedures reduce the run times by at least one order 
of magnitude. For RSSP with time, fuel, and risk constraints, as well as two altitudes, 
our B/B algorithm solves problem instances with 10 by 10 cells and a time horizon 
of 40 typically within 20 minutes. 
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We develop a B/B algorithm and two heuristics (the static bound heuristic 
and the cross-entropy (CE) heuristic) to solve MSP. Among these algorithms, the CE 
heuristic performs better for a broad range of problem instances. 

We also focus on MHSP and develop an exact outer approximation (OA) 
algorithms using several novel cutting planes (multi-cut, strong-cut, relative-cut, and 
symmetric-cut). In empirical studies, this algorithm provides search plans that are 
guaranteed to be within 5% of an optimal plan in less than about 20 minutes for 
problem instances involving 15 by 15 cells, 3 searchers, and a time horizon of 16. 
Instances with 5 searchers are solved even faster. In addition, we prove that under 
certain assumptions the nonlinear multiple homogenous searcher problem is equivalent 
to a large-scale linear mixed-integer program. 

B. FUTURE RESEARCH 

The CE heuristic appears to quickly generate a good feasible solution for MSP. 
However, it does not guarantee solution quality. A possible future area of research 
would be to develop a hybrid algorithm that uses the CE heuristic to generate an 
initial feasible solution for an OA algorithm. Currently we implement the CE heuristic 
using C-|--|- and execute the OA algorithms using GAMS with the CPLEX solver. The 
development of an integrated computational program is worthy to effectively solve the 
MSP/MHSP. Eurthermore, an improved implementation of the OA algorithms outside 
GAMS may reduce cut generation time and prove worthwhile. 

This dissertation does not consider the multiple searcher problem with resource 
constraints (MRSP). In practical applications, the mission planner would like to hnd 
an optimal search plan for multiple “resource-constrained” searchers. A combined 
approach of the RSSP algorithm (Lagrangian directional static bound) with the MSP 
algorithm (CE heuristic) could be developed for this purpose. Eor example, such 
an approach is applicable for the case where the multiple searchers start the search 
operation all at the same time and hnish their mission at the same time. 
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